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ABSTRACT 


An analytical method for predicting the extent of damage in ship collisions 1s developed. 
The method calculates both the longitudinal and transverse extents of damage for the 
struck ship as a function of collision scenario parameters. such as ship speeds. relative 
courses, and collision impact point. as well as considering structural details of the ship. 


This prediction method, or collision model, is used in a Monte Carlo analysis with 
probability density functions defining the specific collision scenario parameters for each 
“case. The Monte Carlo analvsis generates a statistically significant number of collision 
events and results which are applied directly to calculate oil outflow, and from which 
resultant pdf's for longitudinal and transverse extent of damage are calculated to compare 
structural concepts. 


The collision model scenario inputs are initially "calibrated" using a MARPOL single- 
hull model by minimizing the difference between the model result damage pdf s and the 
pdf s specified by the IMO. A quantitative comparison is made between structural 
models for an intermediate oil-tight deck (or "mid-deck ) tanker, and a series of double- 
hull tankers based on calculated oil outflow parameters. 


Of the three ship designs studied. the double-hull series shows the best performance. 
followed closely by the mid-deck tanker. The single-hull ship results predict both more 
frequent and larger spills. 
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1. Introduction 


1.1.Motivation 


The majority of the world’s petroleum travels by sea before reaching the ultimate 
consumer. While this oil is seaborne, it presents a risk to the marine ecosystem and 
shoreline communities. Occasionally. this risk 1s highlighted by spectacular and highly 
publicized events such as the grounding of Exxon Valdez. More often, it is demonstrated 


through less sensational incidents which result in smaller releases of oil. 


In an effort to mitigate this risk, the International Maritime Organization developed 
regulations and associated guidelines which use a probabilistic method to evaluate the 
performance of tankships in groundings and collisions. Alternatives to double-hull ship 
designs must demonstrate performance at least equivalent to double hull "reference" 
designs included in the regulation. This regulation provides an objective method of 


assessing the suitability of potential ship designs, but it has some weaknesses: 


e The current regulation does not consider structural detail. The onlv design attribute 
important to the analysis 1s subdivision. It 1s clear that some structural detail bevond 
this 1s important in determining the response of the structure in these accidents. An 
improved regulation would account for differences in structural detail important to 


structural response. 


e The current method of analvsis assumes that the transverse and longitudinal extents of 
damage are independent. This ts an inappropriate assumption in the case of collision. 
An improved regulation would eliminate this assumption or provide a method to 


incorporate the coupling between longitudinal and transverse damage extents. 


e The current regulation is based on collision data from single hull ships. The current 


regulation uses probability density functions for damage extent to calculate values of 
extreme and mean outflow, and the probability of zero outflow. These pdf s were 


developed from data collected from actual ship collisions, but the data onlv includes 


jl 


collisions involving single hull ships. Applying a pdf developed from single-hull 
ships to a ship with substantially different characteristics may unfairly penalize the 
new ship design by not accounting for improved crash resistance. An improved 
regulation would use pdf's that are developed specifically for the design types under 


consideration, or eliminate the use of damage extent pdf's altogether. 


e The current regulation onlv includes cases where the hull envelope was breached. 


Undoubtedly, there have been instances of collision where no hull failure has 
occurred, and including these cases in the damage extent pdf's may have important 
effects on the probability of zero outflow and value of mean outflow. An improved 
regulation would include the effects of cases where the collision did not result in hull 


failure. 


e The current regulation uses pdf s generated from a small sample size. The pdf's in 


the current regulation were developed from only 52 collisions. Creating a probability 
distribution from such a limited data set is an arbitrary process because there are 
many distribution functions that could fit the data equally well. The final pdf s 
chosen may not accurately reflect the actual probabilistic distributions. An improved 


regulation would base pdf s on a larger sample size, or not use pdf s at all. 


e The current regulation assumes bigger ships suffer more extensive damage. Because 


of the way the regulation uses pdf s normalized bv the struck ship’s length and beam, 
the resulting damage extents are larger when applied to bigger ships. The extent of 
damage should be related to the energy absorbed in the collision, not the size of the 


ship suffering damage. 


To improve the current regulation, a method must be developed which can rationally 
predict the performance of a specific ship design in a given collision scenario. Once this 
is done, the method should be used to predict the performance of the ship design in a 
large number of scenarios. The specifics of each scenario should be determined 
randomly. Analysis of a large number of collision scenarios provides pdf's for extent of 


damage and oil outflow that directly incorporate both the probability that a given 
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collision will result in no outflow, and the coupled nature of longitudinal and transverse 
damage. This also properly credits designs with enhanced resistance to collision. This 
thesis develops such a method, and applies it to several ship designs to demonstrate how 


the current regulations might be improved. 


In a larger view, the proper assessment of risk requires a total system approach [7] where 
the risk is evaluated by separately assessing both the likelihood of an adverse event (such 
as collision, grounding. accidental operational discharge, etc.) and the consequence of 
such events. This thesis attempts to help fill in a missing piece of the overall risk picture 
by providing a rational and quantitative method of assessing the probable collision 
damage to a ship and resulting oil outflow, given the appropriate input distributions 
derived from the specific waterway characteristics. This application of the method is not 


shown here, but would be a straightforward application of the model developed. 





1.2.Review of Collision Analysis Methods 


Previous studies of this subject have taken four different general approaches: finite 
element analyses, model tests, analyses based on the “first principles” of structural 
mechanics, and empirical studies of actual collisions. Each approach has strengths and 
weaknesses that make it suitable for some applications and unsuitable for others. Each 


approach is reviewed briefly to assess suitability for the use at hand. 


Finite element analvses use sophisticated computer programs to construct models of the 
ships involved in a collision. Bv dividing the model structure into many small elements, 
each of which has behavior dictated by its material properties, the structural response to 
given loads can be calculated. The entire collision process can be modeled with high 
precision, including as much structural detail as desired. The drawback to these methods 
Is that development of the finite element models is time consuming, and evaluating the 


collision process is computationally intensive. 


“Model” tests consist of the construction of scale or full size models, usually of onlv part 
of a ship's structure. A collision test is conducted using either a wedge-shaped object or 
another model of a ship's bow as the striking piece. These tests allow carefully 
controlled initial conditions and extensive post-collision analvsis. It is also possible to 
instrument the model hulls with strain gauges or accelerometers to collect time series data 
during the collision process. The drawbacks to this method are that only one collision 
scenario can be studied per model. and that constructing models is both time consuming 
and expensive compared to analvtical methods. Also. if more economical scale models 
are used. issues related to scaling effects arise. These effects range from the scale effects 


of added mass to the scale effects of grain size in the material used for model 


construction. 


Progress is being made toward calculating collision effects from the “first principles” of 
mechanics. Examples of this approach include works by McDermott. et al. [15], 
Reckling [16], and most recently by Wierzbicki [6]. These approaches emphasize closed 


form analytical solutions to a variety of sub-problems within the global perspective of 


- 
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ship collision. For example, solutions have been proposed for plate tearing, plate 
fracture, and failure and resistance of transverse members. These methods and solutions 
are promising, and predict deformation/energy absorption characteristics well in 
laboratory tests. They are difficult to apply in the analysis of a “real-world” collision 
because of the complex interrelationships between modes of failure of each structural 


member. 


V.U. Minorsky conducted the first and best known of the empirical collision studies [2]. 
Since then others have re-validated Minorsky”s original analysis [3] and extended it to 
the analysis of other ship types [4]. The method consists of relating the energy dissipated 
in a collision event to the volume of damaged structure. Actual collisions in which the 
ship speeds, collision angles, and extent of damage are known are used to empirically 
determine a proportionality constant. This constant relates damage volume to energy 
dissipation. In the original analysis the collision 1s assumed to be totally inelastic, and 
motion is limited to a single degree of freedom. Under these assumptions. a closed form 
solution for damaged volume can be obtained, and using the known structural details of 


the ships, extent of damage can be calculated. 


Table 1 summarizes the strengths and weaknesses of these four methods. 


Method Advantages Disadvantages Suitable? 


Finite Element | Precision Complex 
Time consuming model 
construction 

Computationally intensive 
Accuracy strongly dependent on 
appropriateness of modeling 
Expensive 
One collision per model 
Scaling problems 


First Principles | Generality of approach | Still in development 
Basis in physical law 
Empirical Simplicity Oversimplifies collision process 


Table 1: Summary of collision analysis methods 










































Control of conditions 
Extensive analysis 


Model Tests 
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A thorough review of current research and methodologies was conducted by the Ship 
Structures Committee [1]. After evaluating each of the analysis methods described 
above, they concluded that the most promising possibility was to extend Minorsky’s 
original analysis of high-energy collisions by including consideration of shell membrane 


energy absorption. This is the approach taken here. 


a 





1.3. Overview of this Analysis 


This section presents an overview of the major elements of this analysis. Each is 


described in more detail in later sections. 


1.3.1. Collision Analysis Method 


Of the four collision analvsis methods discussed above, only one is suitable for use here. 
The objectives of this work are to develop an analysis method that includes the effect of 
important structural detail, yet 1s simple enough to be readily applied to a large number of 
collision cases and a number of different ship designs ina relatively short time, so that a 


statistically valid distribution of results is obtained. 


The finite element method of analysis is unsuitable for this application because it requires 
a substantial amount of time to develop the finite element model. Finite element analvses 
are computationally intensive, and require a substantial amount of time to complete the 
calculations for even a single collision. Developing a statistically valid distribution 
requires thousands of cases. Conducting thousands of finite element studies for each 
proposed ship design is not a practical approach. Also, the intent of the IMO regulations 
is to provide an objective means of approving ship designs before they are built. Finite 
element models could require a potential shipbuilder to expend substantial effort in 
design development (in order to get the level of detail needed for a finite element model) 


only to find that the design 1s not adequate. 


Experimental analvses are critical in the development and validation of analvtic models, 
but are not suitable for this application. The time and effort involved in constructing 
hundreds or thousands of scale models, all correct in structural detail and scale, then 


subjecting them to collision tests, 1s prohibitive. 


The “first principles” approach is ideal for this application, and is uscd in the analysis of 
grounding events with excellent results [8]. The grounding analysis is facilitated by the 


existence of a well-developed computer program that performs the grounding damage 
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assessment for a given set of inputs. Unfortunately, the application of this kind of 
analvsis to the structural elements of ship side shells, bilges, bottoms and inner bottoms 
has only just begun[6], and no such computer program currently exists. Development of 


such a program is beyond the scope of this work. 


The empirical analvsis method 1s best suited for this application for now. The empirical 
method is simple enough to adapt to different ship designs quickly, but is sensitive to 
major structural details. Calculation of results for a single collision can be done rapidlv. 
Calculation of manv collision cases can be done in a reasonable amount of time and 
support development of a statistically valid sample population. This 1s the method 
utilized in this analysis. In particular, the empirical method of Minorsky is used, but the 
method is extended to provide better prediction of damage in low-energy collisions, and 


is generalized to allow for three degrees of freedom of motion for each ship. 


1.3.2. Scenario inputs 


Each collision scenario 1s defined by a set of parameters that establishes the initial 
conditions of the collision event: the collision model then determines the results of the 
collision. The parameters defining the initial conditions are called "scenario inputs. 
They are determined by randomly selecting each input using a pdf that describes the 


range and variability of that parameter. The scenario inputs are: 
e speed (independently chosen for each ship) 

e collision angle 

e bow entrance angle (for the striking ship) 

e Minorsky energy coefficient 


e initial collision contact point 


striking ship mass 


The pdf' s describing the range and variation for these inputs are developed through a 


combination of data collected from actual collisions and expert opinion, and then 
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modified by a calibration process (described more fully in Section 3.1) to produce 
reasonable results. Figure 1 shows how these inputs describe the developing collision 
scenario. The struck ship’s course angle (wl in Figure 1) is zero at the start of every 


collision. 





Figure 1: Collision scenario definitions 


1.3.3. Collision Kinematics and Simulation 


Minorsky's empirical analysis [2] is generalized to allow for freedom of motion in the x- 
y plane, and rotation about the z-axis for each ship. One unfortunate result of this 
generalization is that a closed-form analytical solution is no longer obtained. Minorsky 
assumed that the only energy important to the collision was that from striking ship 
motion perpendicular to the struck ship. He also assumed that the collision was totally 
inelastic, so that the ships remained joined after the collision. In the generalized model 
all of the kinetic energy is considered and the ships are both free to rotate. This allows 
energy to go into rotational motion as well as deformation of struck ship structure. These 
additional “unknowns” prevent solution of a set of simultaneous equations, and require a 


zi 


time-domain collision simulation to calculate the final state variables (linear and 
rotational velocity) and damage extents. The collision simulation is further discussed in 


Section III. 


1.3.4. Calibration of Input pdf's 


A “calibration” process is conducted to complete the definition of scenario inputs. The 
concept underlying the calibration process is that 1f the collision simulation model is 
accurate, use of scenario inputs drawn from pdf's that describe the “real world” range and 
variation of these inputs should produce result pdf’s that match “real world” damage data. 
Although several shortcomings of the IMO regulation pdfs have been identified, they 
were developed from data collected from actual collisions, and the raw data that forms 
the basis of the pdf's is also available [12]. This makes them a good source of “real 
world" statistics for comparison. Two of the previously described limitations — 
considering only single-hull ship collisions and only considering cases where the outer 
shell is ruptured — are compensated for by using similar restrictions in the calibration 
process. The input pdf's are calibrated to minimize the difference between the damage 
extent pdfs produced by the model and the IMO guideline pdf's. This process is 


described in detail in Section VI. 


1.3.5. Calculation of Damage Extent and Oil Outflow pdf's 


After the scenario input pdf s are calibrated against the current regulatory damage extent 
pdf's, the collision simulation model is run using several dE. ship designs. These 
ship designs were developed using an American Bureau of Shipping (ABS) — proprietary 
software product called “SafeHull”. Each design was verified to meet ABS requirements 
for section modulus. The ships include a MARPOL 73/78 single-hull design (used in the 
calibration process), a tanker with an intermediate oil-tight deck (or *mid-deck" design), 
and several variations on a double-hull design. The individual ship designs are described 


in Section IV. 
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The collision simulation model produces pdf s describing oil outflow and extent of 
damage for each of the ship designs. The oil outflow for each case is calculated directly 
in the simulation rather than from the damage pdf's. The oil outflow pdf's are used to 
determine the probability of zero outflow and mean outflow for each ship type. The pdf s 
describing extent of damage are used only to assess the coupling between transverse and 
longitudinal extents of damage. and to compare the crashworthiness of the different 
designs. The following process description and Figure 2 show the overall approach taken 


for this investigation, and how the calibration process fits into the overall evaluation 


scheme. 
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Overview of the Selected Method 


Choose a means of modeling the structural response. 
Make initial estimates for input parameter pdf's. 
Evaluate a number of collision events for a MARPOL single-hull ship. 


Extract cases resulting in breach of the outer hull, and compare the resulting damage 


extent pdfs to those found in the IMO regulations. 


Modify the input parameter pdf s to provide the best match to the given IMO pdf s 


and underlying data. 
Repeat steps 3 - 6 until the match is satisfactory, or no further improvement is noted. 
Fix the input pdf s. 


Evaluate a statistically significant number of collision events for the ship design types 


of interest. 1.e., single-hull, double-hull and mid-deck tankers. 


Analyze the results. 


Develop structural 






analysis method 
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2. Collision Model Development 


This section describes the time-domain simulation, and the assumptions and 


approximations used to model the force mechanisms. 


2.1. Requirements and Generalization of the Minorsky Method 


Previous work in the area of this studv recommended that Minorsky's analysis be 
extended to include the effects of hull envelope on energy absorbtion [1], and outlined a 
method to generalize Minorsky s analysis to allow for ship motion in more than one 
direction [4]. Both of these concepts are incorporated in this collision model. The 


considerations that guided the development of the model are allowances for: 


forward (surge) motion of each ship 


e lateral (sway) motion of each ship 

e rotation of each ship about it’s own vertical axis (yaw) 

e energy absorption of hull membrane prior to fracture or rupture 

e energy absorption of interior longitudinal bulkheads prior to fracture or rupture 
e calculation of oil outflow resulting from each collision 


e calculation of the longitudinal and transverse extent of the damaged region at each 


moment of time during the collision process 


The first three considerations result in the need to calculate and include the effects of 
hydrodynamic added mass. The fourth and fifth considerations require a model of energy 
absorption for shell plating and bulkheads. The last two considerations require that the 


final collision process model be capable of tracking the extent of damage throughout the 


collision, thereby producing a final damage plan. Each of these requirements 1s 


addressed. 


2.2.Assumptions 


Implicit and explicit assumptions include: 


1. The presence ofa free surface is neglected. No wave interactions or effect of free 


surface on hydrodynamic behavior 1s modeled. 


2. The striking ship's bow is modeled as a rigid triangular structure. None of the 
collision energy is absorbed by the striking ship's structure. This assumption is 


conservative, in that it causes the model to overestimate damage to the struck ship. 


The collision 1s assumed to occur with little to no warning. so no collision avoidance 


UJ) 


actions are taken. This 1s implicitly included in the model by 
e setting initial yaw rate for both ships to zero. and 
e neglecting any effects of backing down on ship velocities 


Assuming that no collision avoidance actions are taken prior to or during the collision is 


also conservative. 


4. The analysis is not a survivability calculation. No calculation of hull girder residual 
strength or damaged stability 1s conducted. The assumption is that the struck ship 


remains intact and continues to float. 


5. Two force mechanisms are considered. One, herein called the "Minorsky force", 
results from plastic deformation of structure internal to the struck ship. The other. 
referred to as the “membrane force’, results from energy absorbed bv plastic 
deformation of the struck ship’s hull or internal longitudinal bulkheads prior to 


fracture or rupture. Each of these forces has an associated direction of action and 


reaction. 
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e The Minorsky force acts in a direction parallel to the direction of relative motion 
between the ships at the point of impact. 1.e., in a direction to oppose relative 


motion. 


e The membrane force acts in a direction perpendicular to the axis of the struck 


ship, including any rotation that occurs during the collision process. 


2.3.Calculation of Added Mass 


As the ship moves through the water, it interacts with the fluid medium. The effect of the 
surrounding water must be accounted for when calculating the result of external forces 
applied to the ship. This is done by including a hydrodynamic “added mass”. Added 
mass is a tensor, the components of which depend on the geometrv of the ship hull, and 


the direction of movement with respect to the ship axes. The form of the tensor 1s: 


where the subscripts indicate the principal body-fixed axes of the ship. The subscript “1” 
corresponds to the principal longitudinal axis and motions in surge. The subscript “2” 
corresponds to the transverse axis and motions in sway, and the subscript “3” to the 
vertical axis and motions in heave. In the absence of detailed information about the hull 
forms and appendages, it is assumed that the ships are symmetric across all three planes. 
This eliminates coupling between the cross-terms in the added mass tensor, which can 


then be written: 


an 0 0 
A= 0 az Q 
0 0 a 


Finally. since this analysis does not consider motion in the z-axis direction, none of the 


elements of the tensor relating to heave motion are retained. This leaves the final form: 


The added mass in surge is calculated by approximating the added mass as equal to that 


of a flat plate having the same area as the midship section of the ship. The added mass of 
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a flat plate for motion normal to the plate is equal to the displaced mass of the 


€ircuniscribed circle [op aenererore: 
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where the midship section coefficient is assumed to be equal to one. For the ships 
examined in this study. this formulation results in an added mass in surge of 


approximately 9% of the ship's displacement. 


A strip theory approach is used to calculate added mass in sway. The coefficient of 
added mass for a two-dimensional rectangular section with side length of 2a under lateral 


acceleration 1s [5]: 
mx) = 4.754 pa" 


Neglecting any effects due to the free surface, and setting the calculated draft equal to 2a 


yields: 
max) =1.189 pT” 


Using strip theory to obtain the added mass coefficient for the three-dimensional ship: 
L 
Uo [max 
0 


OT 


1.189 pT?L 
A 


Co = 


This results in an added mass in sway of about 42.5% of the ship displacement. This 1s 
likely greater than what would be measured for a real ship. This is because the derivation 
treats the ship as a long rectangular structure, neglects the reduction in area at the bow 
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and stern, and also neglects the effect of the free surface. As a comparison, Minorsky's 
original analysis assumed a value of 40% [2]. Experiments have measured this quantity 
at 40% for short collisions [17], but noted that the value is strongly dependent on the 
collision duration. The formulation outlined above is retained because there is no better 
means of calculating added mass in sway without more detailed knowledge of the hull 
and appendage geon:etry. A more accurate approach requires substantially more 
information about hull geometry than exists for the particular ship designs under 
consideration. It would probably not appreciably change the results of the model. Using 
the larger value for the added mass coefficient 1s also conservative because the struck 
ship acquires less translational velocity during the collision, and therefore more energy 


goes into structural deformation. 


2.4. Rotational Inertia 


In addition to linear translation (in the x- and v- axis directions), rotational degrees of 
freedom must be considered for each ship. This requires a calculation of both physical 
mass moment of inertia and added mass moment of inertia about the z-axis for each ship. 


The physical mass moment of inertia 1s calculated by the familiar formula: 
Iss = Tr dm(r) 


where the integration is performed over the entire mass of the body. The notation 
"dm(r)" indicates that the mass distribution is a function of the distance "r^ from the 
origin. For this calculation. the mass of the ship is assumed to be evenly distributed 
along the length of the ship. This is also a conservative estimate, in that the actual mass 
distribution will have relatively less mass at the extreme ends of the ship. and therefore 
have a lower moment of inertia. The model will therefore have less energy going into 


rotational motion, and more into structural deformation. 


The added mass moment of inertia is calculated in a similar manner, using the two- 


dimensional added mass in sway for each incremental section about the axis of rotation: 


2 
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fea = m dmx(r) 


The virtual mass moment of inertia for the ship is then the sum of the physical and added 


mass moments of inertia. 


loor = [eo + feo 


The virtual mass moment of inertia is used to calculate the rotational accelerations and 


velocities resulting from the forces developed during the collision. 


2.5. Overall Energy Balance, and Elapsed Time during collision 


In order to provide an additional check on the final results of the collision process model. 
an overall energy balance 1s developed. By making simplifying assumptions similar to 
those made by Minorsky. an upper bound is calculated for the energy absorbed in 


structural deformation. The specific assumptions for this energy balance are: 


e each collision is totally inelastic, so that after the collision the ships translate as a 


single body, rotating about a common center of mass 
e neither ship has any initial rotational velocity 


e the final physical arrangement of the two ships is such that the striking ship is 
embedded to a depth of one-half the beam into the struck ship, with an angle between 
the two ship’s longitudinal axes equal to the initial collision angle. This sets the final 


mass distribution so the virtual mass moment of inertia can be calculated for the two- 


ship system. 


Making these assumptions, and using conservation of linear and rotational momentum 


yields an equation for the energy absorbed by the ship structure: 
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where 


My = the virtual mass (physical plus added) of ship # 

Va = the initial velocity of ship 7 

Ve - the final velocity of the two ship system 

Gi the virtual mass moment of inertial of the two ship system 
Of = the final rotational velocity of the two ship system 


This energy is compared to the energy absorbed by the structure in the collision model. 
The energy calculated in the energy balance represents an upper bound on the energy that 
could be absorbed by the structure, because the energy balance permits less energy to go 


into other degrees of freedom. 


AS previously mentioned, the collision process model used is a time domain simulation of 
the collision process. In order to execute such a simulation, an appropriate time step 

must be selected. Hutchison [4] analyzed the effects of time step size on solution 
accuracy for a similar time domain collision analvsis. The solution accuracy was judged 
by how well energy was conserved through the collision process by comparing the sum 
of energy absorbed and final kinetic energy to the initial kinetic energv of the system. 

The results of that analysis indicated that the accuracy of the simulation converged 
rapidly as time step size was reduced to a value of T/400, where T is the total time that 
elapses during the collision process. Rather than conducting a similar analvsis for this 
simulation process, a smaller value of T/1000 is used. T is calculated as demonstrated by 


Hutchison: 
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where 
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the aggregate structural thickness, in inches 
the striking ship bow half-entrance angle 


the virtual mass (physical plus added) of ship # 


Again. the energy balance equation is used only as a check on the results of the collision 


process model, and plays no role in the process model itself. The calculation of collision 


duration is onlv used to help select an appropriate time step for the simulation process. 


This shortens the computation time for a given collision scenario by matching the time 


step with the duration of that collision. An equally valid but less elegant approach is to 


choose a conservatively small time step and use that for all collision scenarios. 


2.6.Energy Absorption 


Development of this collision model is guided by previous work done by the Ship 
Structures Committee [1]. They categorized the relative importance of possible energy 
dissipation mechanisms as “primary (or significant), secondary (or not very important), 


: E AT : : 
and tertiary (or negligible)" . The various mechanisms ranked by category were: 


Primary: 


1. Membrane tension in plating, deck and stiffeners 


N 


Plastic bending in plating, deck and stiffeners 


UJ 


Plastic energy in shearing deformation of web frames 
4. Rigid body motion (translation and rotation) 
Secondary: 

5. Elastic bending 

6. Elastic vibration 

Tertiarv: 


7. Thermal 


In this collision process model the secondary and tertiary mechanisms are assumed to be 
negligible. The primary mechanisms are combined into two categories. The first 
category describes the energy absorbed by membrane tension in hull and internal 


longitudinal bulkhead plating. This 1s referred to as “membrane energy along witha 


' The quote comes directly from Volume 1 of Reference [1]. The original phrasing is 
kept to convey the panel's intent to portray even secondary mechanisms as “not very 


important’. 


corresponding “membrane force”. The second category covers all other primary 
structural interactions and deformations. This is accomplished by relating the amount of 
energy absorbed to a volume of structure damaged via an empirically determined 


coefficient. This is referred to as “Minorsky energy” and the “Minorsky force”. 


2.6.1. Membrane Energy 


Minorsky's original analysis of ship collisions forms the basis of the collision process 
model developed here, but Minorsky concentrated on “high-energy” — collisions where 
the ships involved are large and have high initial relative velocities. This focus produced 
a good linear relationship between volume of damaged structure and absorbed energy for 


high-energy collisions, but did not produce good results for low energy collisions. 
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Figure 2: Original Minorsky Correlation” 


: Adapted from reference [2]. 


Noting the area of poorly fit data points in the low energy regime of 


Figure 2, the Ship Structures Committee recommended that Minorsky’s method be 
extended into the low-energy regime by adding some consideration of the energy 
absorption capacity of the hull envelope [1]. Reardon and Sprung accomplish this by 
using a theoretical model of a wedge cutting a plate [5, 13, and 14]. The revalidated 
relationship shows almost exactly the same numerical value for the Minorsky coefficient, 


but improves the fit in the low energy regime. See Figure 3. 


Revalidated Minorsky Relation 
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Figure 3: Revalidated and extended Minorsky correlation 


In Appendix C of Reference [1]. a different method to do this 1s outlined by Jones, and 
further elaborated by Van Mater. This method treats the panels of the hull envelope as 
thin, broad “beams”, with “pinned-pinned” end boundary conditions. The ends of the 


beam correspond to the panel attachments to internal ship framing. 


> Adapted from reference [3]. 


Used together with Reardon and Sprung’s results, this method allows: 
ə calculation of the energy used in deforming the beam 

e arbitrary location of the application force between frames 

e an approximation of the deflection prior to plate fracture or rupture 


For this study it 1s also important to consider any inner hull membranes (as in the case of 
double-hull and mid-deck tankers) and internal longitudinal bulkheads. This is important 
not onlv from an energy absorption and crashworthiness point of view, but also when 

determining whether or not a given cargo bulkhead has ruptured and allowed oil outflow. 
Because of these issues, the extension proposed by Jones and Van Mater is applied to all 


of the vertical longitudinal plate structures. 
The specific formulas used in these calculations are from [1]: 


atBw? a 
Eai = M99 ev 0M 
L b 
Wlinit — 0.452a 


where 


ee energy absorbed by the deformed membrane 
Oy = yield stress of steel 
t thickness of plate 


B = effective breadth of plate 
W = deflection of plate 
L = length of plate between clamped ends 


a&b = shorter and longer distance from point of force application to 
clamped ends 
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Figure 4: Membrane force mechanism concept 


Figure 4 shows how the deflection of the plate is calculated from geometrical 


considerations. 


The plate breadth (B) 1s not the vertical extent of the plate, but an effective breadth. This 
effective breadth is calculated by requiring the expectation value of energy absorption to 


match the shell energy absorption value of 28.4 MJ calculated by Reardon and Sprung in 


[3]: 


The expectation value is calculated over the possible range of the short leg length, a. 
This provides an effective breadth that will result in a expectation value for Emem of 28.4 


MJ, as calculated in [3]. 


In each time step, the plate deflection 1s calculated from the geometry established by the 


bow shape and ship positions. Then there are two possibilities: 


If the deflection is less than the deflection limit (wig) given above, the energy absorbed 
by the plate is calculated. Dividing this work done in the time step by the incremental 


deflection in the time step gives the magnitude of the force due to the plate deflection. 


4] 


The direction of the force is assumed to be perpendicular to the longitudinal axis of the 


struck ship, and forms an action/reaction force pair, which acts on both ships. 


If the deflection limit 1s exceeded, the energy absorbed by the membrane is set to zero, 
and the calculated force is therefore also zero. This is also used to define the rupture 

status of this membrane. In the simulation code, this triggers a transition to a different 
subroutine that no longer calculates plate deflection or energy absorption until and if a 


new plate or membrane structure is encountered. 


2.6.2. Minorsky Energy 


The other structural forces are combined into a single mechanism based on Minorsky’s 


relationship between energy absorption and volume of damaged structure. 


Figure 5: Minorsky mechanism concept 


This is implemented by calculating the volume of newly damaged structure at each time 
step. Since one of the initial assumptions 1s that the striking ship s bow is a rigid triangle, 


this becomes a relatively simple problem in trigonometry. See Figure 5. 


To calculate the newly damaged structure, the area covered by the triangular bow’s 
incursion into the struck ship s side 1s determined, and then the area covered during the 
previous time step is subtracted. Finally, the newly damaged area due to longitudinal 
relative motion is calculated and added to give the newly damaged area during each time 
step. This area is multiplied by the Minorsky coefficient to give the energy absorbed 
during the time step. This is converted to a force by dividing the energy by the distance of 
relative travel during the time step. The force is modeled as an action/reaction force pair 


directed to oppose the relative motion. 


3. Collision Scenarios 


3.1. Selection and Calibration of Input pdf's 


The collision simulation model is a deterministic process model, i.e., if the simulation is 
run repeatedly with the same initial input parameters, the same output will be produced. 
In order to explore the probabilistic nature of the damage extents resulting from collision, 
a stochastic process 1s employed to choose the parameters that define the initial state 
variables of the system. These parameters that define the initial conditions of an 
individual collision scenario are called "input parameters". They are drawn randomly 
from a set of probability density functions, each of which describes the relative likelihood 


of occurrence of any particular range of values for that parameter. 


3.2.Ship Speeds 


The ship speeds were selected from a pdf defined by a bimodal distribution made up of 
two normal distributions with mean values of 5 and 10 knots, each with variance of 1 


knot?. 





Target PDF 
+ + Calculated Distribution 


0.25 
02 

= 
un 
= 
v 
O 
> 0.15, 
5 
oO 
de: 
e 
Ses 

0.05 

0 M 
2 4 6 8 10 12 14 
Ship Speed (kt) 


Figure 6: Ship speed histogram and pdf 
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This distribution was selected to match that used by Rawson [8], and is based on rational 
argument and expert opinion. The assumption is that tankers spend the majority of their 


time operating in one of two modes: 
e Transit mode — represented by the 10 knot mode of the distribution 
e Maneuvering mode — represented by the 5 knot mode of the distribution 


The variance of each mode is an estimate based on the opinion of the author and several 
professional mariners. Figure 6 1s a histogram of collision speeds selected for one of the 


ships (based on 5000 collision cases) compared to the speed probability density function. 


3.3. Collision Angle 


Collision angle is defined as the angle of incidence between the ships at the moment of 
impact. This input parameter 1s based on a uniform distribution. The choice to use this 
distribution is based on rational argument. Different distributions could be postulated (c.f. 
[9]) but would depend on a particular route and waterway. As an example, a route that 
includes a long straight channel would result in a collision angle pdf with higher density 


near zero degrees and 180 degrees. 
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Figure 7: Collision angle probability density function 
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A route that includes an extended port approach with many turns and blind spots would 
result in a collision angle pdf with higher density near 90 degrees. The proposed uniform 
pdf represents a compromise for a generic tanker in worldwide trade. Collisions 

occurring at a relative angle of zero degrees are constrained to have an initial impact 
point at the bow of the struck ship. Collisions occurring at a relative angle of 180 degrees 
are constrained to have an impact point at the stern, and are only allowed if the striking 
ship's speed exceeds the struck vessel's speed. Figure 7 shows a histogram of selected 
collision angles (based on 5000 collisions) compared to the specified pdf from which 


thev are drawn. 


3.4. Impact Point 


The point along the struck ship where the striking ship's bow initially makes contact is 
called the impact point. This point is allowed to vary with equal probability along the 
entire length ofthe ship. The selection of this pdf is based on rational argument. The pdf 
for impact point is simply a linear function with density equal to one along the length of 
the ship. Figure 8 shows a histogram of impact points based on 5000 collision cases 


compared to the pdf from which the impact points were drawn. 
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Figure 8: Initial impact point probability density function 
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The impact point includes parts of the ship forward of the forward-most cargo bulkhead, 
and aft of the aft-most cargo bulkhead. The collision dynamics account for the effect of 
forces and moments applied at these locations. The extent of damage is recorded and 
compared to the actual cargo tank boundaries in order to calculate oil outflow. This 


maximizes the “realism” of the collision simulation, rather than constraining the collision 


to begin within the cargo block. 


3.5.Minorsky Coefficient 


For each collision scenario generated, a particular value is selected for use in the 
Minorskv relationship between energy absorption and volume of structure damaged. The 
pdf used to select this value is a normal distribution with mean equal to 47.1 MJ/m3 and 
standard deviation of 8.8 MJ/m3. This distribution is based on a validation of Minorskv's 
original work done by Reardon and Sprung [3], including the addition of new data points 
from collisions that have occurred since Minorsky s work in 1959. Figure 9 shows a 5000 
case histogram of this parameter, along with the pdf from which the selections were 


drawn. 
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Figure 9: Minorsky constant pdf 
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3.6. Bow Entrance Angle 


The shape of the bow of the striking ship 1s important because it determines the volume 


of structure subject to damage during the collision. 
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Figure 10: Bow half-entrance angle probability density function 


In this analysis, the shape of the striking ship's bow is idealized as a triangle, with no 
rake. For each collision scenario, the bow half-entrance angle is selected using a pdf with 
a normal distribution, with a mean value of 38 degrees and standard deviation of 5 
degrees. This distribution 1s based on data presented in [3] and [9] showing 
representative bow half-entrance angles for a range of ship displacements, and 


adjustments made during pdf calibration. 


Figure 10 shows a histogram of selected bow half-entrance angles (based on 5000 


collisions) compared to the specified pdf from which they are drawn. 


3.7. Striking Ship Displacement 


For each collision scenario, the particulars describing the struck ship are given and 
constant. The striking ship is not known, and its characteristics must be chosen using a 
distribution as for other scenario input parameters. A common approach to this problem 
is to assume that the striking ship and the struck ship are identical in all respects. This is 
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based on the assumption that “like ships” travel the same waterways (being engaged in 
the same trade), and are therefore more likely to have collisions with other similar ships. 
This approach is not satisfying here, because the amount of energy the striking ship 
imparts to the collision process (and therefore the extent of damage) is strongly 
dependent on the mass of the striking vessel. The mass of the striking vessel should be 
chosen from a distribution that, like collision angle, reflects the waterway environment in 
which the ship is or will be operating. For this study. which is not specific to any 
particular waterway, the striking ship's displacement was selected from a normal 
distribution with a mean of 150,000 metric tons, and a standard deviation of 30,000 
metric tons. This choice for this distribution was based on data from [9]. and validated 


by the calibration process. The distribution 1s shown in Figure 11. 
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Figure 11: Striking ship mass histogram and pdf 


3.8. Calibration of Input Parameter pdf's 


To ensure that the results of the simulation are reasonable, the damage extent pdf s 
resulting from the collision simulation model and scenario inputs were compared to the 
equivalent pdf s in Regulation 15F. The pdf's in the Regulation are based on data 


collected from real collisions involving single hull vessels during the period 1980 — 1990. 
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They only include collisions in which the hull envelope was breached. In order to base 
this comparison on similar data sets, the collision simulation model was run using a 
MARPOL single hull ship as the struck vessel, and the output data was discarded for 
collisions where the outer hull is not ruptured. This initial comparison was favorable, and 
provided a basis for confidence in the validity of the model. Input pdf's required only 
minor adjustment. The end result of the calibration process was a set of input parameter 
pdf s that are verified to give reasonable results in this model, and may be used by others 


conducting similar analvses. 


Following the initial comparison. a calibration process was undertaken to improve the 
correlation between the pdf s in the Regulations and those produced bv the model. The 
correlation is measured by constructing a “goodness of fit^ parameter similar to the “R- 
squared" parameter used to quantify correlation in a linear regression method (c.f., [11], 
especially Chapter 9). After the collision model was run and inappropriate cases 
discarded, the damage extent pdf's were calculated. The simulation pdf values (over 
discrete ranges) were compared to the IMO pdf over that same range. The difference is 
squared, and the sum taken over all the discrete ranges in the pdf, then divided by the 


number of ranges. 


In equation form: 


> Di, (Sa E Ro) 


n 


R 


where 
Spdf = the value of the simulation pdf over the ith range 
Rpdf = the value of the IMO pdf calculated over the ith range 
n = the number of ranges taken 


and the sum is taken over all n ranges. This value was calculated for each of the pdf's in 
the Regulation: longitudinal extent of damage, extent of transverse penetration, and 
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longitudinal location of the center of the damaged area. An overall goodness-of-fit 


(hereafter “fit”) was calculated by summing all squared differences and dividing by the 


total number of sample ranges. 


Some of the input parameter pdf s were not considered for modification. The 
characteristics of the ship speed pdf's are considered fixed, since they are based on work 
done concurrently by Rawson [8] and provide good results in a grounding analysis. The 
pdf describing initial point of contact for the collision is also fixed, since only a uniform 
random distribution reproduces the regulatory pdf describing the location of the center of 
the damaged area. The parameters left to vary were those describing the striking ship’s 
bow half-entrance angle, and the striking ship displacement. The medians of both these 
distributions were varied over reasonable ranges in an attempt to maximize the fit 


parameter. The results are shown in Figure 12. 
Calibration surface plot 
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Figure 12: Calibration sensitivity plot 


This calibration surface was generated by conducting nine simulation runs with the 
parameters varied over the range of interest. The data was then fit using a cubic 
interpolation scheme. It is apparent from the data collected that the local maximum for 
goodness of fit lies somewhere in the middle of the range explored. The median values 
for these distributions were selected bv this process, and are as described previously. The 


scale used in Figure 12 makes the difference in fit over the explored range seem large. 
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In fact, the best fit parameter obtained was 0.7754, and the worst fit was 0.7411, a 
difference of only about 5%. This means that the results of the collision simulation 
model are relatively insensitive to variations in these two parameters within a reasonable 


range. 


3.9. Single Hull Results Compared to IMO pdf's 


3.9.1. Center of Damaged Extent 


The pdf describing the location of the center of the damaged section corresponds fairly 

well with the pdf specified in the Regulation. There is some random variation around the 
uniform density level of one. Superimposed on this random variation is a bias toward the 
stern of the ship. This bias 1s caused by the motion of the stricken ship, which is assumed 


to be in the forward direction. 
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Figure 13: Location of center of damaged section for single hull tanker 


3.9.2. Longitudinal Extent of Damage 


The longitudinal extent of damage predicted by the collision simulation model matches 


the Regulation very closely. One difference is that the model predicts some damage that 
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exceeds 0.3L, which is the upper limit of the Regulation pdf. The frequency of these 
cases is low. The transverse extent of damage associated with these cases is also low 
(this can be seen from Figure 23.) These damage patterns result from high energy 


collisions with low angles of incidence. 
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Figure 14: Longitudinal Extent of Damage for Single hull Tanker 


Another difference is that extremely short longitudinal extents do not appear as 


frequently as the IMO Regulation would predict. 


3.9.3. Transverse Extent of Damage 


The pdf describing transverse extent of damage differs substantially from the regulation 
pdf. The most notable characteristic is the collection of damage cases with transverse 
extent around 0.2B. The collision simulation model produces this result because of the 
longitudinal bulkhead located at 0.1875B. The additional resistance presented by this 


bulkhead stops many collisions, and reduces the relative velocity so that many more are 


Di 
G2 


stopped in the next few meters. The model also predicts a small probability of transverse 


damage exceeding 0.3B, whereas the Regulation never predicts this. 
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Figure 15: Transverse Extent of Damage for Single hull Tanker 


This particular pdf causes the most difficulty in the calibration process. No variation of 
input parameters will eliminate this "spike" from the model output. This is the limiting 


factor in improving correspondence with the IMO pdf's. 
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4. Ship Designs 


4.1.General Specifications 


A family of representative tankers was designed by Rawson [10] for calibrating the input 
scenario pdf's and estimating the effect of structural enhancements on crashworthiness. 
The tankers include a MARPOL single hull tanker. five double hull tanker variants and 
an intermediate oil-tight deck (mid-deck) tanker, all of Suezmax (150,000 dwt) 
dimensions. The single hull tanker 1s designed consistent in material and configuration 
with vessels in service between 1980 and 1990, the period included in the data compiled 
by the classification societies to generate the current IMO damage pdf's. This design is 
used to calibrate the scenario probability densitv functions by matching the calculated 
damage extent density functions to the densitv functions provided in the Guidelines. The 
double hull and mid-deck configurations are designed using current shipbuilding 


practices and used for comparisons between design alternatives. 


All designs have the same pri:cipal dimensions listed in Table 2, with bulkheads located 
to maintain equal cargo capacities and compliance with MARPOL Regulations for 
protective location of segregated ballast tanks. maximum tank volumes and double hull 
requirements. Scantlings are the minimum allowed by current classification societies 
standards, as determined by the American Bureau of Shipping's SafeHull system. The 
effect of structural enhancements on crashworthiness 1s studied using five separate 
double hull variants. Each variant is a derivative of the original baseline double hull 
model. with either the plating thickness. stiffener sizes, stiffener spacing or frame spacing 
modified. For each new variant design, the remaining structural parameters are re- 
examined using SafeHull to ensure optimum (1.e., minimum) compliance with 


classification scantling requirements. 


Table 2 and Figure 16 provide the general specifications to which each ship configuration 
is designed. Details are provided in Figure 17 through Figure 19 and Table 5 through 
Table 10. In each figure, a plan view and transverse section are shown. Scantlings are 


listed in the tables. The plate thickness listed is the average of the plate thickness over 
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the breadth of the respective bottom, side, deck or bulkhead. Frame spacing is +.7 m for 


all ships except the final double-hull variant, where frame spacing is 3.7 m. 


Length Between Perpendiculars 
Beam 

Draft 

Depth 

Double-Bottom Depth 

Wing Ballast Tank Width 
Displacement 

Deadweight Tonnage 

Plating material 

Cargo Tank Arrangement 


Single-Hull 
264 m 

48 m 

16.8 m 

24 m 

N/A 

Om 

178,411 mton 
~ 150k 

MS24 

5X3 

plus 2 slop tanks 


Double-Hull 
264 m 

48 m 

16.8 m 

24 m 

2.4 m 

Zim 

178,411 mton 
~ 150k 

HT36 

6x3 

plus 2 slop tanks 


Table 2: General Specifications of Ship Designs 


Figure 16: Typical tankship profile 


IOTD 

264 m 

48 m 

16.8 m 

24 m 

N/A 

515515011 

178,411 mton 
~ 150k 

Lo 

6 X 2 OVer 6 al 
plus 2 slop tanks 


4.2. Single Hull 
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Figure 17: Single Hull Plan and Section 


Slop Tank Cargo #5  Cargo#4  Cargo#3 Cargo#2 Cargo #l 
or Ballast or Ballast or Ballast 
Port 2195 | 6520 8878 A E 5975 
enter N/A 20057 DOSE cM 29041 19918 
Starboard 2195 6520 8878 8878 8712 3013 


Table 3: Cargo Tank volumes for Single-hull ship design (m3) 


The scantling dimensions used 1n the Minorsky and membrane force calculations for the 


single-hull ship are: 





Component Thickness 











Side shell plate thickness 25170 Cm 
Bottom plate thickness 2 lelem 
Upper Deck plate thickness 2:2 Cm 

Aggregate deck thickness 4.381 cm 
Internal longitudinal bulkhead plate 1.952 cm 


thickness 


Table 4: Single Hull Scantlings for collision analysis 


4.3. Double Hull Ships 
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Figure 18: Double hull Plan and Section 


Slop Tank Cargo #6 . Cargo Z5 Cargo#4 Cargo #3 


or Ballast 
Port T955 8767 8767 8767 8767 
Center N/A 16908 13828 13828 13828 
Starboard 19553 8767 8767 8767 8767 


Cargo #2 
or Ballast 
8348 
13167 
8548 


Table 5: Cargo Tank volumes for Double-hull ship designs (m3) 


4.3.1. DH - Baseline 


Cargo #1 


O 00 Un 
UU») Un UG 
NO — NO 
NO ON. NO 


Rawson developed five separate versions of the double-hull configuration. The baseline 


DH design represents an “optimized” ship that uses the minimum weight of steel to meet 


ABS requirements. The scantling dimensions used in the Minorsky and membrane force 


calculations for the baseline double-hull (DH) ship are listed in Table 6: 


Component Thickness 
Side shell plate thickness 1.8cm 
Bottom plate thickness 1.881 cm 
Inner Bottom plate thickness I: eei 
Upper Deck plate thickness 2 dcm 
Aggregate deck thickness S1 0/512. (6106. 
Inner skin plate thickness 1.842 em 
Internal longitudinal bulkhead plate leo) Cina 


thickness 


Table 6: DH Scantlings for collision analysis 


4.3.2. DH1 


DH variant #1 (DHI) is derived from the baseline by increasing plate thickness to 150% 


of its original value. The 


Component Thickness 


Side shell plate thickness ecm 
Bottom plate thickness 2.822 cm 
Inner Bottom plate thickness 2.656 cm 
Upper Deck plate thickness 3.0 cm 
Aggregate deck thickness 8.478 cm 
Inner skin plate thickness 207 CHI 
Internal longitudinal bulkhead plate 2:359 emt 
thickness 


Table 7: DHI scantlings for collision analysis 


dor ELI 


DH variant 22 (DH2) is derived from the baseline double-hull by increasing the 
scantlings of all stiffeners so that their contribution to total section modulus is 150% of 
the original design. Because the collision model used here does not consider the impact 
of individual stiffeners, and the plate thickness 1s not changed from the original design, 
there is no difference in collision performance between the baseline ship and DH2 in this 


model. This variant is not discussed further. 


4.3.4. DH3 
DH variant 23 (DH3) 1s derived from the original double-hull design by reducing the 


stiffener spacing to 75% of that used in the original design. 


60 





Thickness 
1.614 cm 
1.405 cm 





Component 
Side shell plate thickness 
Bottom plate thickness 









Inner Bottom plate thickness 1.4 cm 

Upper Deck plate thickness ST 
Aggregate deck thickness 4.66 cm 
Inner skin plate thickness 1.694 cm 
Internal longitudinal bulkhead plate 1.256 cm 


thickness 





Table 6: DH» scantlings for collision analysis 


Plate thickness is reduced to the maximum extent possible while still meeting ABS 


requirements for section modulus. The dimensions used in the collision analysis for DH3 


are shown in Table 8. 


4.3.5. DH4 


DH variant £4 uses the same scantlings as the baseline ship, but the frame spacing is 
reduced from 4.7 m to 3.7 m. In the collision simulation model frame spacing is used to 


determine the dimensions of plating for calculating membrane energy and force. 


61 


4.4. Intermediate Oil-Tight Deck Ship 






Siop Tk —s 







Slop Tk — 






Figure 19: Mid Deck Plan and Section 


Slop Tank Cargo #6 Cargo#5 Cargo#4/ Cargo #3/ Cargo #2/ Cargo #l 


o Ballast Ballast Ballast 
5 cep ee ee 
Starboard 1626 7526 7326 7326 7320 6909 3979 
Lower 17903 14632 14652 14652 13819 7959 


Table 9: Cargo Tank volumes for IOTD ship design (m3) 


The scantling dimensions used in the Minorsky and membrane force calculations for the 


IOTD ship are: 


Component Thickness 


Side shell plate thickness 1.8 cm 
Bottom plate thickness 1.8 cm 
Inner Bottom plate thickness 1.56 cm 
Upper Deck plate thickness pM 
Aggregate deck thickness 5.68 cm 
Inner Skin plate thickness 1.845 cm 


Centerline Bulkhead plate thickness — 1.622 cm 


Table 10: Mid Deck scantlings for collision analysis 


The mid-deck tanker represents an alternative design under the Regulation. This ship has 
shell scantlings similar to the MARPOL single hull tanker, but the decks and bottom are 
reduced because of the presence of the internal horizontal mid-deck. This ship also has a 
centerline bulkhead which extends from the upper deck to the mid-deck, but not to the 
inner bottom. Another distinctive feature of this design is the rather wide ballast tankage 
outboard (double sides). This functions as protectively located ballast and provides a 


measure of protection in collisions. 





5. Results 


5.1.Mean Outflow and Probability of Zero Outflow 


Mean oil outflow and the probability of zero outflow are calculated for each ship design. 


The results are shown in Table 11. 


The probability of zero outflow calculated in this analysis is not the same as the 
probability of zero outflow as used in the Regulation. In the Regulation the calculated 
value is a conditional probability that is more properly described as “the probability of 
zero outflow given a collision that results in hull rupture.” The value calculated in this 
study is also a conditional probability. It is "the probability of zero outflow given a 


collision’. 


The value of mean outflow is non-dimensionalized by dividing the outflow volume by 


the total cargo volume of the ship. 


Ship Design Probabilitv of Zero Outflow 


Single Hull 
Mid-Deck 
Double Hull (baseline) 


Double Hull (enhanced plate) 
Double Hull (reduced stiffener 
spacing) | 
Double Hull (reduced frame spacing) 





Table 11: Probability of zero outflow and mean outflow results 


Retemimne te Dable di: 


The mid-deck design has the highest probability of zero outflow, but also the highest 


mean outflow. The high mean outflow is related to the subdivision scheme chosen by the 
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designer. Once the cargo boundary is breached, 75% of the oil in that cargo section is 
lost. This ship would have substantially improved performance if an intermediate oil- 
tight deck were combined with a more typical "three tank across” arrangement or a center 
line bulkhead in the lower tank. These arrangements would also reduce potential intact 
stability problems associated with free surface effects during loading and unloading. To 
explore this, the simulation was run again with a lower centerline bulkhead added to the 


design. Pp remained constant at 62%,but mean outflow dropped from 10% to 8%. 


The double hull designs all show roughly similar performance. The design with 
enhanced plate thickness shows the best performance within the double-hull group. This 
is expected from the way the Minorsky method relies on in-plane elements for energy 
absorption. The double hull designs have a lower Po in collision, due to the relatively 
small protective layer of the double-side, but have the lowest mean outflow because of 


their greater subdivision 


The single hull ship shows Po second only to the mid-deck. and mean outflow between 
mid-deck and double hull. The mean outflow of the single hull is also adversely 


impacted by the chosen subdivision. 


5.2. Extent of Damage pdf's 


The extent of damage pdf s presented are not conditional on rupture of the ship S muii 
Because of the difference in conditionality. these pdf s should not be compared to the 


IMO pdf's, only between ship designs analvzed here. 


5.2.1. Longitudinal Extent of Damage 


The pdf describing longitudinal center of damage 1s roughly uniform, as expected. There 
is a slight bias toward the stern of the ship. This ts a result of the forward motion of the 
struck ship. The initial point of impact was selected from a uniform pdf. This sets one 
end of the damaged extent. The location of the other end depends on the relative speeds 


of the two ships. The only way the “end” of the damaged section can be forward of the 
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initial collision point is for the striking ship to have more velocity in the x-direction that 


the struck ship. This tends to shift the center of damage aft. 


+ Calculated Distribution 
12 Uniform Distribution 


eo 
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Longitudinal center of damaged section (Loc/LBP) 


Figure 20: Probability density function for center of damage 


The effect is not dramatic, but it is seen consistently through all of the simulations. A 
plot of the pdf from the Regulat ons and the pdf resulting from the simulation for the 
single hull ship is shown in Figure 20. The results for the other ships are 
indistinguishable from the single hull case. Since the oil outflow is calculated from the 
damage plan for each individual collision, this pdf does not play any part in assessing the 
performance of the ship designs. Pdf's for the other ship designs are not presented 


because they are essentially identical. 


9.2.2. Longitudinal Extent of Damage 


The pdf s describing longitudinal extents of damage are shown in Figure 21. 
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Figure 21: Longitudinal extent of damage pdf's 


These pdf's are similar within the group because all the ships have similar side protection 


svstems. The double hull actually has the least protection in side collision. The single 


hull has substantial protectively located ballast tankage, and the mid-deck has even more 


side protection. 


Note that longitudinal damage extents exceeding 0.3L are predicted for all ship designs. 


The Regulation does not predict any damage beyond this length. There are three 


explanations for this: 


The small number of collision cases that form the basis of the Regulation pdf's did 
not show damages of this extent. It is possible that a larger sample size would have 
shown this. The pdf s developed in this analysis are based on 5000 collisions, and the 


pdfs show that damage exceeding 0.3L 1s uncommon. 
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e The assumptions that are included in the simulation model are intentionally 
conservative. This produces greater extents of damage than a model without these 
assumptions. The conservative assumptions that cause this include the treatment of 


added mass and the rigid striking bow assumption. 


e The Regulation pdf's are based on collisions that result in rupture of the hull. Some 
of the collisions that result in damage extents exceeding 0.3L do not result in hull 
rupture. They can be characterized as “glancing blow” collisions. Collisions of this 


type are excluded from the IMO pdf's. 


Overall, correspondence with the “real world” data is good, and shows that the collision 


simulation model produces reasonable results. 


5.2.3. Transverse Extent of Damage 


The pdf s describing transverse extent of damage for each ship design are shown in 


Figure 22. 
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Figure 22: Transverse penetration pdf's 


The general trend shows two main concentrations of transverse penetration. There is a 
cluster of collisions that are halted or nearly halted bv the shell ofthe ship. There is 
another cluster of cases that are stopped by the next internal membrane. This can be seen 


in all ship designs. 


e The single hull ship shows this second cluster around 0.2B, which corresponds to the 


longitudinal bulkhead at 0.1875B. 


e The double hull ships show this second cluster around 0.3B, which corresponds to the 
internal longitudinal bulkhead at 0.2975B. It is not possible to see a similar grouping 
between the inner and outer hull because they are so closely spaced. Two meters is a 


difference of 0.042B. This is too fine to be resolved by the bin sizes of 0.05B. 
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e The mid-deck ship has an internal cargo bulkhead at 0.1145B so the second cluster is 
easier to identify. The mid-deck ship also has another cluster at 0.5B, which 


corresponds to the centerline bulkhead. 


The mid-deck and double hull ships generally show greater transverse penetration. This 
is because they have less in-plane structure (deck thickness) to absorb energy via the 
Minorsky mechanism. All the ships were designed to satisfy ABS requirements for 
section modulus, but the extra material used in constructing the double hull or internal 


oil-tight deck allowed thinner plate to be used in the bottoms and decks. 


DH1 shows transverse penetration comparable, but still greater than the single hull ship. 
This is seen in spite of having greater aggregate deck thickness (5.572 cm vs 4.381 cm 
for the single hull). This shows that shell plating is a significant energy absorber. The 
significant energy absorption by shell plating 1s what causes the data scatter in the low 


energy regime of Minorsky’s original analysis. 


Thinking about collision resistance in terms of the Minorsky interaction, it is clear that 
the more efficiently a ship is designed (assuming traditional structural designs). the less 
collision resistance the ship will have. This is similar to the argument that improvements 
in engineering knowledge have resulted in decreased “safety room” as the design margins 


have been whittled away over time by improved knowledge of structural response. 


5.2.4. Joint Longitudinal/Transverse Damage pdf's 


Figure 23 through Figure 28 are joint probability density functions showing the 
distribution of coupled transverse and longitudinal extents of damage. The damage 
results of each collision scenario are recorded and analyzed as a set so that the transverse 
and longitudinal extents can be plotted as a dependent pair. Each ship shows slightly 
different results, mostly in transverse extent of damage. The differences result from the 
same factors previously discussed relating to transverse extent of damage. It is clear that 
there is coupling between transverse and longitudinal extent of damage. One can see 


how the pdf for one extent (say transverse extent) depends on the other extent of damage 


D 


by imagining a two-dimensional “slice” of the joint pdf taken at a representative value of 
the other extent (longitudinal, in this case). As one moves the “slice” along the selected 
axis, the shape of the slice changes. It is impossible to capture this effect with the 
methods in the current Regulation without a joint pdf plot. The method used in this 
analysis captures the effect completely. This coupling effect may or may not be 
important. The only way to tell 1s to use the outflow characteristics developed by this 
analysis to calculate a “pollution prevention index” for each of these ships, and then 


compare to the pollution prevention index calculated via the IMO Regulation. 


Joint PDF for Transverse and Longitudinal Extents of Damage 
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Figure 23: Joint damage pdf for single hull 


The single hull joint pdf shows the impact of the internal longitudinal bulkhead at 
0.1875B. Asa ship strikes, the first resistance encountered is from the shell. In this 
model, the direction of the force developed by the shell tends to “turn away” the 
incoming ship, and directs collision damage into longitudinal extent. If a ship "punches 


through" the exterior shell, damage tends to be transverse until a second longitudinal 


membrane is contacted. The same turning force is exerted, and damage is directed 


longitudinally. 


Joint PDF for Transverse and Longitudinal Extents of Damage 
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Figure 24: Joint damage pdf for mid-deck tanker 


Figure 24 shows the same effects, including the protective bulkhead at 0.1142B. This is 
seen as a smaller ridge inboard and parallel to the ridge created by the shell membrane. 
Of the cases where a striking ship proceeds past this bulkhead, nany proceed all the way 
to the center line bulkhead. This ship design has the most transverse penetration of those 


analvzed. 
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Figure 25: Joint damage pdf for baseline double hull 


Figure 25 through Figure 28 show the joint damage pdf's for the double hull series. They 
have the same characteristics because of their similarity in design. They show similar 
behavior to the other ships with respect to the longitudinal membranes. The double hull 
series does not show the effect of the inner hull very clearly because of the grid size used 
in the joint probability plots. This is a trade-off. A fine mesh size would provide the 
resolution needed to see the effect of the closely spaced hulls. A larger mesh size 
prevents random variation in the bin populations from obscuring the trends. The mesh 
size used here is the smallest that prevents random vanation from becoming problematic. 


Using a finer mesh requires more than 5000 cases. 
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Figure 26: Joint pdf for double hull DH1 
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Figure 27: Joint damage pdf for double hull DH3 
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Figure 28: Joint damage pdf for double hull DH4 
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5.3. Oi] Outflow pdf's 


The oil outflow pdf's for all ships are shown in Figure 29 through Figure 34. 
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Figure 29: Oil outflow pdf for single hull 
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Figure 30: Oil outflow pdf for mid-deck 
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Figure 31: Oil outflow pdf for baseline double hull 
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Figure 32: Oil outflow pdf for DH] 
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Figure 33: Oil outflow pdf for DH3 
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Figure 34: Oil outflow pdf for DH4 
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Examination of these pdf's reveals three items of interest: 


e The first column on the left-hand side of the graph corresponds to the probability of 
zero outflow. The way the pdf graphs are constructed, this probability can be 
approximated from the value of the “zero” column, although some cases of very small 


outflow are included in this column. 


e The height of the remaining columns corresponds to the likelihood of an oil outflow 


of a given size. 


e The discrete nature of the oil outflow shows that there are only a few values of oil 
outflow possible. The particular values depend on the ship design, and are a result of 


assuming that once a cargo tank is ruptured all the oil in that tank is lost. 


These figures are interesting because they depict the full range of possible results given 
that a collision has occurred. The probability of zero outflow includes cases where the 
hull is not ruptured as well as those cases where the hull fails. The remainder of the 

graph shows the range and probability of all outflows. The current regulation does not 


require such a plot, although one could be constructed using the Regulation methodology. 
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6. Conclusions 


A rational method of calculating the probabilities of oil outflow has been demonstrated. 


The method: 


Considers structural detail 


Treats ships appropriately based on size 


Considers the coupled nature of longitudinal and transverse extents of damage 


e Is tailored to the particular ship design 
e Hasa statistically significant basis for prediction 


This method is rapidly adaptable to other ship designs. Itis fast, simple, and treats 
structural differences in a rational manner. It is an improvement over the current IMO 


Regulation methodology. 


With regard to the particular ship designs considered here, the following conclusions are 


drawn: 


The double-hull provides the best performance of the designs considered 


e The mid-deck design is superior to the double-hull in terms of providing maximum 


chance of preventing all outflow in a collision, but has higher mean outflow. 


e The single hull design shows larger spills than the double-hull and more frequent oil 


spills than the mid-deck. 


e Subdivision is critical in limiting oil outflow. The mid-deck design could be 
comparable to the double-hull if an improved subdivision scheme were implemented. 
The performance of the single hull ship would also be improved if more subdivision 


1 


were used. 
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7. Future work 


There are several areas where this analysis could be improved. Most of the 


improvements lie in the elimination or refinement of assumptions made in this work. 


Efforts that would provide the greatest benefit include: 


Elimination of the rigid bow assumption. Damage surveys following collisions 
typically reveal that a substantial amount of damage is also done to the bow of the 
striking ship. This analysis has assumed that the striking bow is impervious to 
damage. This causes more energy to be absorbed by the struck ship. and therefore 


produces more extensive damage. 


Include resistance mechanisms for major transverse elements. This analysis does not 
explicitly consider the transverse elements. Finding a way to include the effects of 
transverse webs and bulkheads would improve the responsiveness of the model to 
structural detail, and eliminate another conservative assumption that tends to over- 


predict extents of damage. 


Collect and implement better statistics for input parameters. The pdf's for input 
parameters used here are the result of expert opinion and the calibration process 
described in Section VI.. This process was consciously undertaken with the goal of 
producing damage extent pdf's similar to those in the Regulation. A better approach 
would be to collect enough data from ship operations to base the input pdf s on real- 
world statistical data. A comparison should then be made to real-world statistical 


data on actual collisions. The collision data should not be “scrubbed” to include only 


the collisions that result in hull failure. The collision data should be sorted by ship 


type. 


~ 


The model should be run using actual collision data to see how well it predicts the 
results of a single case where the inputs and results are well known. This would 


guide further improvement efforts. 


The model’s use of the Minorsky mechanism should be eliminated. This mechanism 
represents a range of failure modes aggregated into a single constant. When first- 
principles methods are mature enough to represent the majority of thes: failure 
modes, the individual failure mode calculations should replace the broad-brush 


Minorsky method. 
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9. Appendix 


All of the computer code used in this analysis is included. All of the code is in the form 


of MATLABO script files. MATLABO Version 3 was used throughout this analysis 


The script files are presented in the order used. The first script, calculate.m, calls the 


following scripts as required to complete the calculations and produce the output. 


To run the code, all files need to be in a common directory in the MATLAB path. 
Typing “calculate” at the MATLAB prompt starts the execution. The analysis is guided 


by user input prompted by a series of questions. 


In order to analyze different double-hull ships, changes need to be made in the 
“calculate.m” script for new deck thicknesses, shell plating thickness, and internal 


bulkhead thickness. 


In order to change the frame spacing, a change needs to be made in “constants.m” in the 


variable “spacing”. This is marked with comments in the code. 
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% calculate.m 


Ene Contre! Script to run eacgssubrouteue 4 035025 
Date created: 10/15/97 


o? o9 oe 


Last revision: 3/10/98 
% Inputs: Nümber oft desired. suns. ne 
OUtpuc: 


o? oe 


Results of momentum and energy balance for final 
velocities, collision energy, and time step 


% "dt" - an estimate of the appropriate time step to be 
used in the time domain simulation 

eit 

clear 

Ide 

n - input('How many runs, please? E 

Syoe = inputt'Press 1 for single hull, 2 foredousue-uj eco 
LOST ET 


% Set up storage arrays 

E =zeros(n,1); % pre-allocated memory for EA 

Eshell = zeros(n,1); *tpre-allocated memory for Emem 
Vl=zeros(n,1); % pre-allocated memory for Vell 
V2=zeros(n,1); % pre-allocated memory for Vel2 


bow alpha = zerosin, 1); + presallocarted mens euis 
Len - zeros(n,1); * pre-allocated memory for L 

R = zeros(n,1);% pre-allocated memory for R - a flag for 
outer shell rupture 

Cor =" Zeros ay) ; Y pre-allocated memory for collision 
angle (phi) 

P = zeros(n,1);% pre-allocated memory for penetration depth 
Cen = zeros (n, l); % pre-allocated memory for center of 
damaged section 

ICP = zeros(n,1); % pre-allocated memory for initial 
contact point 

FCP = zeros(n,1); % pre-allocated memory for final contact 
pole 

Minorsky = zeros(n,1); % pre-allocated memory for the 
Minorsky constant 

Mass = zeros(n,1); % preallocated memory for striking ship 
mass 

Time = zeros(n,1); % pre-allocated memory to track collision 
time 

outflow = zeros (1); % pre-allocated memory for oil 
outflow 

test = 0; 

Lif type == 1 
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test = input('Press 1 if this is SoD e e SE 
ASA E ae 
end 
constants 
EZ Eros (o 
ona ms 


vargen 
energy 
% if Ea <0 
2 Ea = 0; 
% end 
LL CU 
SBH = [13,54.1,95.2,136.3,177.4,204.7,218.5,265]; 
singlel 
if ruptUre o 00 
single2 
single3 
if rupt úre? po 
single4 
end 
end 


oil = (Tank1P + Tank1S) * 15935 + Tank2 * 29041 + 
(Tank3 + Tank4) * 29592 4 Tank5P * 6520 4 Tank5 * 290529 
Slop iis ee HS 
oil = oil/166577; 
end 
1 ME SI 
sigma y = 3.6e+08; % yield stress of HT36 steam 
in Pa 
KE = GKE* (36/24); % correction Lor mi 
vice MS24 per Daidola & Pet 


t s Omer % aggregate deck thickness for DH1 
t plate = .01614; % shell thickness for DH1 
DBE = [137 46779,112 7145, 178) 201 gue > eee 
Coup wed 
i eee O 
t plate- TOEA * inner skin thickness 
doubl z 
double3 
if rupturo 2p 20 
t plate = .01256; % centerline pukki PEE 
thickness 
double4 
doubles 
lp trot MN NE 
doubles 
end 
end 
ena 


oil = (Tank1P + Tank18) * 9658 + Tank2P * Sedo 
Tank2C * 13167 + (Tank3P + Tank5P + Tank6P) * 8767 + 
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a oc rank 4. Tank5C) * 13828 s Tank6Cc *—$69935 
Soot UOS. 
oil = 011/164079; 


end 
if type == 3 
sigma y = 3.6e+08;  *yreld stress of HT36 steel, 
in Pa 
GKE = dKE* (36/24); % correction for HT36 steel 
vice MS24 per Daidola & Pet 
tum 05687 * sets aggregate deck 
thickness for the IOTD design 
Pao O * shell thickness for 
TOTD 
EBR € [13746;79,112,145,178 21129749 EDS F 
Woje oli 
LAO ES > 00 
e placen mOoi % inner skin thickness, 
in meters 
rocd: 
iotd3 
Sp eU E UD E 0 
E plater Toe 7 MESA e 
bulkhead thickness, in meters 
iotd4 
rOEdS 
=> o O 
LOCE 
end 
end 
end 
oil = (Tank1P + Tank1S) * 3979 + Tank1BOT * 7959 + 


(Tank2P + Tank2S) * 6909 + Tank2BOT * 13819 + (Tank3P + 
Tank3S + Tank4P + Tank4S + Tank5P + Tank5S F Tanki 
Tank6S) * 7326 + (Tank3BOT + Tank4BOT + Tank5BOT) * 14652 + 
Tank6BOT * 17903 + SlopTk * 1626; 
odl qd 16279 

end 
write 
end 
* Remove cases that fail the integration error test from the 
data set 


howmany 
% zero out the cases that fail the test 
tar e in 
qu Time (IIASA) 
DA = 0; 
Emus 
ees O 
bow alpha(i) = 0; 
DA RO 
He Em Or 
end 
if P(i) > Beam 
ues ese 


EI 


end 


end 


$ Now, remove those cases from the population using the 
"honzemesd rurc/e en 


Len = nonzeros (Len) ; 

P = nonzeros(P); 

Cen = nonzeros (Cen) ; 

bow alpha = nonzeros (bow alpha) ; 
CA = nonzeros(CA) ; 

ICP = nonzeros (ICP): 


% Outflow is not removed here because applying the 


Nem Ze nese 
% function would also take out all the cases where the run 


Was 


6 acceptable but no outflow occurred. These zeros will not 
affect 
% the calculated mean as long as the sum is divided by the 
proper 


5$ population size. 


* Non-dimensionalize the output vectors 
$Len - Len/LBP; 

P = P/Beam; 

Cen = Cen/LBP; 


o? 


Plot a summary figure 


figure(1) 
Cite 
bone o ae 
SubplGt (3, 1,1) 
he ewer 
hist (Cen, bin) 
xlabel ('Longitudinal center of damaged section!) 
ylabel('number of occurrences!) 
bin=[0.05:0.1:1]; 
SUBO OS 2) 
Holda on 
hist (Len,bin) 
xlabel ('Longitudinal Extent of damaged section') 
ylabel('number of occurrences!) 
SUPE LOC gem) 
22 


Þin= (0 1:0.05:1]; 

Ie cum 

hist (POIN) 

xlabel ('Transverse Penetration') 
ylabel('number of occurrences') 


*% Produce other output plots, and goodness-of-fit data if 
chis 2S a calibration run 
Ares 
calibrate 
else 
QMieDuUE 
end 


SO LO ESAS pdf 
point 
rotate3d 


7o Vargen.m 


5 Routine to generate random variates for use in the 
collision script rour mie 


* Input: the number, n, ot variates to produce. Gurren, 
entered in script as 1000 


Output: Generated variates 


* Variable label Variable description 
% Ul Struck ship speed (kt) 
% U2 Striking ship speed mie) 
% alpha Striking ship half-entrance 
angle (degrees) 
% phi Collision angle (degrees 
relative, from struck ship) 
% L Location of collision point aem 
struck ship (meters from FP) 
% Lnd The collision poine, M eE 
dimensionalized by LBP 
% dKE The value for the energy 
absorption coeffiecient in Minorsky's equation 

m2 The mass of the striking ship 


Date created: 9/12/97 
Last updated: 4/23/98 


o9 o9? o? 


reset the random number generators to new states. 


o? 


rand('state' sumo roOO0*cloek)5 
raundn('state'ssaum lOO Clock, 


% Generate variates in proper distributions and ranges 


% Input for speed generation section 


Dm e: 

lomu = 5; % mean velocity peaks in knots 
hisigma = l; % and standard deviations 
qose = tí: * for high and low speed 


alphamu = 38; striking ship half-entrance angle 


alphasigma = 5; 


if rond >- T 5 


MUSE. 

Sigma = hisigma; 
else 

Mi =. emu 

sigma = losigma; 
end 
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UL =- mu r randn*tsigma; 


rand('stare’,sum(100*clock))>; 
reandmisstate'  sum(i00*clocme) 


me rand > 0.5 


uS DEOS 

sigma - hisigma; 
else 

mu - lomu; 

sigma - losigma; 
end 
U2 = mu + randn*sigma; 


% set impacting bow half-entrance angle 
alpha = alphamu + randn*alphasigma; *%* Normally 
distributed about alphamu 


* set collision angle (phi) 


cheek. =i. 
while check ==1 
Or 6 0 cane 


% phi = 90 + randn*45; % Normally distributed on [0- 
180] degrees 


This version is for a two mode normal dist for phi 


o\? 


Z Nimi = 135; 

% Pomur 5- 

% Same 13 On 

% PE mane >- 2s 

% mu - himu; 

a else 

% mur miu 

% end 

6 pasa aa Sama 


eneek = om (pha i ponhi- i Oor 
end 


* Set collision location point L, in meters 
Lo = LEP*rand; 
ma > TTO LBP, * Lnd is the impact point non- 
dimensionalized by LBP 
qr DNE % Sets collision point for end=on 
COLLIS Ons eyes: 
Lo = LBP; 
end 
Peso hi eel 
LO = 1; 
end 
* select value for the Minorsky resistance function: 
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dkB = 47.1 + ranan*8e.8: 


% Choose a mass for the striking ship in lton, convert to 


kg, and build 
% a physical mass array - this replaces the value assigned 


o 
. 


in Constants. m 


120000 + > Gandn=3 0000. 
mi *o2407 2204 ee 
he Alin or oO 


3 3 
be 
D 


* Calculate the principal dimensions of the striking ship by 
scaling by the cube Geer toe 
% the mass ratio. 


LBP2 = ((m1/m2)^(1/3))*LBP; 
Beam2 = ((m1/m2)~(1/3))*Beam; 
Dral - “(mima sob E 


* By maintaining dimensional similitude, the added mass 
Coet clients rele ln constan m E 
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% energy.m 


% This script calculates the initial energy and momentum of 


the two ships, 
* the final momentum and velocity of the two-ship system 
( 


assumes that they 
5 travel together in the end state), and then the final 


energy of the system. 


* The difference between the initial energy of the system 


and the final energy of the 
% system is used in the next script, which will determine 


how the energy is 
‘6 expended in deformation of the struck ship's structure. 


Date created: 10/7/97 


A A 


Laot rey Salome 4/23/98 
senputse Parameter variables from "vargen.m" 
% Ship masses - from "constants.m" 


Added mass tensor for each ship 


[NEU S Ea - the energy that must be expended in 
structural deformation 
$ dt - an estimate of the appropriate time step 


to be used in the time domain simulation 


* Calculate added mass tensor for each ship 
5$ theta(1 or 2) is the angle between the ship's 
principal axes and the ship's velocity vector. Prien sre 
+ ther collision, this 1S zero ropebset s 


theta 0. 
theta? = 0: 


% Calculate components of added mass tensor for x & y 
Mee akon 


Al = mi*[(alii*sin(thetal*pi/180)^2 + 
a22*Ccos(thetai*piy 160) 2) Mali 
a22)*sin(thetal*pi/180)*cos(thetal*pi/180)); ((all- 
a22)*sin(thetal*pi/180)*cos(thetal*pi/180)), 
(ali*sin(thetal*pi/180)^2 -« a22*cos(thetal*pi/180)^2)]; 

A2 s m2*[(all*sin(theta2*pi/180)^2 + 
a22*cos(theta2*pi/180)^2), ((ali- 
a22)*sin(theta2*pr/I80)*cos(theta2*ps Noc NEM DEDE 
a22)*sin(theta2*pi/180)*cos(theta2*pi/180)), 
(all*sin(theta2*pi/180)^2 « a22*cos(theta2*p1/180)^2)]; 
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* Combine physical and added mass matrices for total 
"virtual" mass matrix 


VM1 M1 + Al; 
VM2 = M2 + A2; 


It 


* Construct9velocity vectors for “esa chest SEE MN E 


meters/sec) 
* so that V1 = (x=velocity, y velocity seme 


jc 1456889009 39) + [01,015 
Vel2 = 1.688* (12/39) * [-U2*cos (phi*pi/180); - 
U2ssin phi pii SO 


* Construct momentum matrices for each ship (units of 


kg*m/sec) 
P1 - VM1*Vell; 
P2 = VM2*Vel2; 


Total momentum is then: 
PT = P1 + P2; 


o? 


oy 


with magnitude 
PF = sare Ei ccc qc es 


o? 


at angle to global coordinate system of: 
chi = men? (PT(2), PT(1)): 
chi deg - chi*180/pi; 


5 With standard "small changes in mass distribution" 
assumption outlined in thesis, the final virtual mass tensor 
iS: 

VMF = VM1 + VM2; 


* and the final translational velocity is then: 
VF = VMF\PT; 


* Rotational Energy Calculation 


% Find center of mass relative to struck ship origin 


I MINOS. 
x = m2*(((LBP2/2) -L) + (LBP/2) *cos (phi*pi/180) ) / (m1+m2) ; 
y = (m2* (LBP/2) *sin(phi*pi/180))/(m1+m2); 


rt ocn M CE EE 
% at angle toe global crirgin or. 


beta = atan2(y,x) ; 
beta_degree = beta*180/pi; 


Calculate physical mass moment of inertia for both ships 
66 le AMES) MZ 
662 = MI (L BP2 2) DF 


o 
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% Calculate added mass moment of inertia for both ships 
(assumes pure sway motion arises 
% from rotation about the final center of mass, so this is 
the same as added mass in sway) 

J66A1 - (2.378*rho*Draft^2*LBP^3)/24; 

J66A2 = (2.378*rho*Draft2”2*LBP2"3)/24; 


% Combine to get virtual mass moment of inertia 
J66Vl1 = J661 + J66A1; 
J66V2 MAA XM OA 


i 


$ Use parallel axis theorem to calculate the virtual mass 
moment of inertia about the 

* System center of mass for the striking ship. Assuming 
again, that the motion of striking ship is 

6 entirely in sway 

% I-new = Icg + Mr“2 


J66V2C 2 J66V2 * m2*(1-«a22)*(((LBP/2)-L-x)^2 + 
MO BR, 2 eos (php Ln aa a 


* Use parallel axis theorem to calculate the virtual mass 
moment of inertia about tne 

5 System center of mass for the struck ship. Assuming, as 
before that the motion of the struck ship is 

entirely in sway 

I-new = Icg « Mr^2 


o? o? 


J66v1lC = J66VI + mi*(l+ag2) = 0 09) 7 C DIE 


* Combine the virtual mass moment of inertia of Ship 1 about 
the system center of mass with the 

* virtual mass moment of inertia of Ship 2 about the system 
center of mass to obtain the 

6 virtual mass moment of inertia of the two-ship system 
about the system center of mass 


JSG O NO goo lc 


* Using this to solve for the final rotational velocity 


[o] 


gives: 


r2 = OGBP/2) - LUE the IEEE Ie 
striking vessel's linear momentum acts 


wf - (1/J66F)*(r2*(sqrt((P2(1)^2 + P2(2)°2)))*sin(phi) 
- rf*PF*sin(beta-chi)); 

wf deg per sec - wf*180/pi; 
* Now, using these quantities, calculate the difference in 
initial and final energy states of 
+ the system. This 1S an approximation er the energy sehen 
must be absorbed by the structure. 
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* It is an approximation because the final virtual mass 
moment of inertia is calculated by 

% assuming that the final mass distribution is the same as 
at moment of impact and under 90 

% degrees collision angle. 


Ea = 1/2*((dot(P1,Vell)) + (dot(P2,Vel2)) - 
(EGE (PT, VE) ) = ODE M VN MC 


* These quantities calculated to assist in troubleshooting 
code - remove for faster execution 

KE1 = 1/2* (dot (P1, Vell) ); 

KE2 1 2HMaot (P2 Velo 

KEE 1/2 COLE VE 

KE TA O SNE WE: 


of oV o9 o? 


% and an approximation of the total time elapsed during the 
colli stems 

T = 
(pi/2)*sqrt((1/(t*300000*tan(alpha*pi/180)))* ((VMEZ(I1,1) 9VMER 
1. T0 2 CUNEINUS 1) 8 2028 NER ESTE 


% and deviating slightly from Hutchison"s work Enee 
stėép for SimulseungNebrsq9 m isrone ds. 


step = 7/200; 


% (Hutchison used T/100) 


100 


% singlel.m 


$ Script to perform time-domain analysis for single-hulled 
tanker collision. 

s The im andicates that this script 1s for phase teervene 
Sollicion. Which 1s 

4 from the time of impact until the shell membrane ruptures. 


N 


censu cC from energy.m 
$ va from energy.m 
% V2 from energy.m 
% VM1 from energy.m 
% VM2 from energy.m 
% alpha from vargern.m 


o 


Output: Generated variates and values 


Date created: 11/3/97 
Last updated: 4/23/98 


o? o 


oe? 


Reset flags for tank breaching to zero: 
Tanke OO, 

Tankis = 0; 

Tank2 
Tank3 = 
Tank4 : 
Tanks 7 
TankSP = 0; 
Sep Tice ee 

%* Determine nearest transverse structures, and distance to 
each for use in applying 

* Van Mater's extension to Jones method.... 


LÀ 
Li 
t 


JA ES 
while BH(j) < L 
D tl 
end 
a= BH(j) - L; 
b BE (ee 
% Ensure variable "a" represents the "short leg" of the 
strained plate 
E a»b 
cea; 
a=b; 
bee 
end 


* Calculate" detlection at which" pUrare duris Xem voee e 


extension to Jones 
eic Mm DS a 
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% Initialize new variables (Subscript 1 = struck ship, 
Subscript 2 = Striking ss ip, 

Eme 205 

X1=0 ; 

Img 0 

X2= (LBP/2) -L+(LBP/2) *cos (phi*pi/180) ; 
Y2= (Beam/2)+(LBP/2)*sin(phi*pi/180); 
omegal=0; 

omega doti qao 

omega2- (phi«180)*(pi/180); 

omega dot2 - 0; 

ce Eo: 

CP = Mel 

T2 <= Vemos 

rupture Too 

rupture2 - 0; 

cuu 

ddepth = 0; 

depth =<08 

Fmin = 0; 

MS O 

Eabs = 0; 

Emem as O 

NOA 


relvel = 1; % a temporary value to get through the first 
cycle of the time-step routine. 
6%%%6%%%S% Begin time-step routine %%%%%%%%% 
while abs(defl) « abs(defl lim) 
if relvel « endvel 
break 
end 
* Calculate new postions and rotations at end of time step 
time = time + step; 
= ee Pia) Sore: 
Yl = Y1 + TI(2)*step, 
X2 = X2 + T2(1)*step; 
Y2 = YO PUO) stem 


omegal = omegal + omega dotl*step; 
omega2 = oméga2 + omega dot2*step; 
[= Lp xam 
Pen = Pen + dPen; 
1f Pen > Pmax 
Pmax = Pen; 
end 
TE Un MB. 
L=LBP; 
break 
end 
Je ee 0) 
[05 
break 
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end 
1f Pen > Beam 
Pen = Beam; 
break 
end 
if Pen < 0 
break 
end 


% Calculate relative translation, penetration, and change in 
impact point in this time step 

S1 = [T1(1)«((LBP/2)-L)*omega dotli*sin(omegal), 
T1(2)-((LBP/2)-L)*omega dotl1*cos(omegal)]; 


S2 IEA = ee 2/2) emegancers*sinvemeqaz ™ 
T2(2)- (LBP2/2) *omega dot2*cos(omega2)] ; 
reltrans = (S2-S1)*step; 


o 


SeOColeulLatemetreckien OL relative translation 
zeta = atan2(reltrans(2),reltrans(1)); 


* Calculate penetration and change in location 

dPen = sqrke ((reltrans( (11) 2 + 
(reltrans(2))^2)*cos(omegal + (3*pi/2) - zeta); 

dL s sqrt((reltrans(1)) 2 + (reltrans(2)) 2) *sin (emegar 
+ (3*pi/2) - zeta); 


* Calculate the membrane deflection 
defl = defl + reltrans(2)+((LBP/2)-L) *sin(omegal) ; 


Calculate the resistance force of the membrane 
Lf derhe<*6 
Emem - sigma y*t plate*breadth*defl1^2/spacing; 
Fmem = Emem/abs (defl); 


o 


else 
break 
end 
Perder Me wder] Lim 
Sa c iim, 
Emem - sigma y*t plate*breadth*defl^2/spacing; 
Fmem = Emem/abs (defl); 
rupe => 
end 
% Calculate the force resulting from the "Minorsky 
mechanism" 
% Rtt is the total "resistance factor" 
% dRt is the differential resistance factor for this 
time step 
% depth is the distance of penetration 
% ddepth is the differential distance of penetration 
* t is the aggregate in-plane structure thickness 
5$ alpha is the bow half-entrance angle 
% This formula accounts for the triangular poes 
Geometry 
% and dynamic collision angle, but places no 
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% limits on striking ship beam. Should be modified so 


that if width exceeds bean, 
* remaining area is rectangular.... 


ddepth = sqrt (reltrans (1) 2 + %seRerans (3) 2 

Rtt = (depth*2) *t*tan(alpha*pi/180) / (1- 
((tan(alpha*pi/180))^2/((tan((omegai-omega2)*pi/180)^2)))); 

dept a = depen Nee c sis 

dRt =  (depth^2) *t*tan(alpha*pi/180)/(1- 
((tan(alpha*pi/180))^2/((tan((omegai-omega2)*pi/180)^2)))) + 
abs(dL)*Pen - Rtt; 


* Calculate the delta-KE from the Minorsky relation; 
Emin - dKE*10^6*dRt; 


* The corresponding force is 
Fmin = Emin/abs (ddepth) ; 


Calculate resulting accelerations from the membrane force 
For Phase I, the membrane force is assumed to act 
perpendicularly to Ene mule si PC. 

of the struck ship, and the Minorsky force acts to oppose 
the direction of relative morton 


o op 


o\? 


Sol (Struek Ship) 
Calculate the virtual mass 
First, calculate the angle of resultant force 


compared to the ship principal axis 


o? o o 


metal) = pi/2; * the membrane force is 
always normalta the Ste as nao X 
% Now, calculate the added mass matrix based on this 
angle 

Al - mi*[(aii*sin(zeta1)^2 + a22*cos(zetal) 2), 
((a11-a22)*sin(zetal)*cos(zetal1)); ((al1i- 
a22)*sin(zetal)*cos(zetal)), (all*sin(zetal) 2 = 
ag2* cos (zeual) 2a 
6 Combine with physical mass to get the virtual mass 
maeri 

VMI = MEF- Al; 
% The acceleration of translation in the global X 
coordinate is: 

aXlmem = Fmem*sin(omegal)/VM1(1,1); 
% The acceleration of translation in the global Y 
coordinate is: 

aYlmem = -Fmem*cos (omegal)/VM1(2,2); 
% Calculate the angular acceleration about the ship c.g.: 

sithe urent contact roin iS: 
CP = [X2 - ((LBP/2)*cos(omega2)), ((LBP/2) - 


L)*sin(omegal)]; 
% so the arm that the force acts through is: 
arm - sqrt((X1-CP(1))^2 « (Yi - CP(2))"2); 
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% so the angular acceleration is: 
omega dotdotimem = -((0.5-Lnd) /abs(0.5- 
Lnd) ) *Fmem*arm/J66V1; 


% Now calculate accelerations due to the Minorsky 


y Nec. 
% The Minorsky force is assumed to act in the direction 


opposite of relative motion. 
“= cuamcecd:hnlcsShoncc Loui amo rrenemt direction we must 


recalculate the virtual mass 


$ First, calculate the angle of resultant force 
compared to the ship principal axis 

zetal - zeta - omega; 5 the angle of 
Minorsky [orce to bue SEES ces 
% Now, calculate the added mass matrix based on this 
angle 

Iu mir (orale tal 2 + a22*cos | zetal 21 
((al1-a22) *sin[(zetal)*cos(zetal)); ((al1- 
a22)*sin(zetal) *cos(zetal)), (all*sin(zetal)“2 + 
el *coszet 2) 
% Combine with physical mass to get the virtual mass 
matrix 


VIRE MI F AL; 
* Now calculate the accelerations due to the Minorsky force: 
aXimin = Fmin*cos (zeta) /VM1(1,1)>;> 
aYlmin = Fmin*sin(zeta) /VM1 (2,2); 
2 dotdotmin-= -| (LBP/2)-L)*Fmintsin(zeťta Tpi 
omegal)/J66V1; 


% Sum the accelerations due to membrane and Minorsky for the 


total accleration due 
fs que to relative motion and interaction Cor Ehiswe we estep 


axl = aXimem + aXlmin; 
aYl = aYlmem + aYlmin; 
omegal dotdot = omega_dotdotimem + omegal dotdotmin; 


Snip. 
Calculate the virtual mass for the membrane force 
acceleration: 


o? o9 


% First, Calculate thewangle of resultant Tord: 
compared to the ship principal axis 

zeta2 = omega2 - (omegal + pi/2) 
% Now, calculate the added mass matrix based on this 
angle 

A2 = m2*[(all*sin(zeta2) "2 + a22*cos(zeta2) 2), 
((aliea22)*siyn(zeta2)*costzeta2) (alle 
a22)*sSsin (zetaZz)*cos(zetal))), Wativsin (ger) 2. 
aZz2*Cos (zetaZ eZ) 
% Combine with physical mass to get the virtual mass 
matrix 


UNAS MZ E AZ 


% The acceleration of translation in the global X 
coordinate is: 
aX2mem = -Fmem*sin(omegal) /VM2(1,1); 
5 The acceleration of translation in the global Y 
coordinate is: 
aY2mem - Fmem*cos (omegal)/VM2(2,2); 
The angular acceleration about the ship c.g. is: 
omega2 dotdotmem = (-Fmem*sin(omegal-omega2- 
1/2) 4 (LB P2720) )/366V2; 


o? 


O 


Now calculate accelerations due to the Minorsky 
interaccion, 

% The Minorsky force is assumed to act in the direction 
Opposite of relative motion. 

% Since this force is in a different direction we must 


recalculate the virtual mass 


o? 


$ First, calculate the angle of resultant force 
compared to the ship principal axis 

zeta2 = pi - omega2 + zeta; % the angle of 
Minorsky force to thefstrik ing onl P 
% Now, calculate the added mass matrix based on this 
angle 

A2 = m2*([(all*sin(zeta2) 2 + a22*cos (zeta2y ae 
((all-a22) *sin(zeta2)*cos(zeta2)); ((al1- 
a22)*sin(zeta2)*cos(zeta2)), (all*sin(zeta2)72 + 
a22*cos(zeta2)^2)]; 
% Combine with physical mass to get the virtual mass 
matrix 


VM2 = M2 + AZ; 

* Now calculate the accelerations due to the Minorsky force: 

aX2min = Fmin*cos(zeta-pi) /VM2(1,1); 

avY2min = Fmin*sin(zeta-p2)y VM2 (2) oe 

omega 2 dorado tal Mee M *Fmin*sin (omega? - zeta) /J66V2 
* Sum the accelerations due to membrane and Minorsky for the 
total accleration due 
* due to relative motion and interaction for this time step 


aX2 = aX2mem + aX2min; 
aY2 = aY2mem + aY2min; 
omega2 dotdot = omega2 dotdotmem + omega2 dotdotmin; 
5 E new velocities 
TECO 1) eer Seer 
A A A E Stee 
T2(1) = T2(1) +%axk2*stee, 
TAI o E 


e omega dotil + comegaldetasissuepr 
omega dot2 = omega dot2 + omega2 T 


telvel = sqrt (veltranms (ree er 
Eabs - Eabs + Emin + Emem - Emem last; 
Emem last = Emem; 

6 Determine which tanks have been breached 
singlecargo 


end 
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% single2.m 


Y Script to perform time-domain analysis for single-hulled 


tanker Collision. 


% The "2" indicates that this script is for phase 2 of the 


Eqllisgteon, walen as 
inner longitudinal cargo 
s bulkhead is contacted. 


o 


* Input: all the dynamic variables from singlel.m 
% at Erom vriteim 

% M tf COM: Wisi m 

% V2 from write.m 

% VM1 from 

% VM2 from 

% alpha from write.m 


* Output: Generated variates and values 


% Date created: 11/3/97 
* Last updated: 4/23/98 


$%%%%%%% Begin time-step routine %%%%%%%%% 
Pen < SS1 

lvel « endvel 

break 


K 
(D (D 


end 
Fabs = Eabs + Emin; 


% from the time the outer shell membrane ruptures until the 


% Calculate new postions and rotations at end of time step 


time = time + step; 
E Xl il) A St ee 
MiS Y ie acter. 
Wu lI IX Tempe 
puc YO LUDERE 
Senegal = Onegai OMedaydort* Seep, 
omega2 = omega2 + omega dot2*step; 
ele dL, E 
Pen = Pen + dPen; 
1f Pen > Pmax 
Pmax = Pen; 
end 
n M 
[n 
break 
end 
if L » LBP 
15 — T1 
break 
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end 
if Pen > Beam 
Pen) = Beam 


break 
end 
if Pen « O 
break 
end 


% Calculate relative translation, penetration, and change in 
impact point in this time step 
% S1 is the total velocity (from lineaxw and 2ctaeweonall 
motion) of the smpact'pount sm M 
+ S2 is the same velocity rem ho ri 

S1 = [T1(1)+((LBP/2) -L) *omega_doti*sin(omegal) , 
T1(2)+((LBP/2) -L) *omega_doti*cos (omegal) ] ; 


S2 = [T2(1) - (LBP2/2)*omega dot2*sin(omega2) , 
T2(2) * (LBP2/2) *omega dot2*cos (omega2)] ; 
reltrans =- 922991055 cu E 


* Calculate direction Otrereleative Lreaas ak 
zeta - atan2(reltrans(2),reltrans(1)); 


* Calculate penetration and change in location 

deren Toart (reltrans (D T2 
(reltrans(2))^2)*cos(omegal + (3*pi/2) - zeta); 

dL - sqrt((reltrans(1))^2 « (reltrans(2)) 2)*sin(omegal 
+ (3401/2) 6t) 


5 Calculate the force resulting from the "Minorsky 
mechanism" 

Rtt is the total "resistance factor" 

dRt is the differential resistance factor for this 
time step 

depth is the distance of penetration 

ddepth is the differential distance of penetration 
t is the aggregate in-plane structure thickness 
alpha is the bow half-entrance angle 

This formula accounts for the triangular bow wedge 


Q9 oV 


o9 o? (T o o o o o 


geometry 
and dynamic collision angle, but places no 
limits on striking ship beam. In the rutcure Acces 


be modified so that if width exceeds beam, 
5$ remaining area is rectangular.... 


ddepth - sqrt (reltrans(1)^2 -« reltrans(2) 2); 

Rtt = (depth^2)*t*tan(alpha*pi/180)/(1i- 
((tan(alpha*pi/180))^2/((tan((omegai-omega2)*pi/180)^2)))); 

depth = depth + ddepth,; 

dRt = (depth*2)*t*tan(alpha*pi/180) /(1- 
((tan(alpha*pi/180))^2/((tan((omegai-omega2)*pi/180)^2))))- 
abs(di)*Pen - Rtt; 
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:xme-dcuatesrhe delta-KE from the Minorsky velerron,; 


o 


Emin = dxkxE*10"6*dRt; 


% The corresponding force is 
Fmin = Emin/abs(ddepth) ; 


c Shy SE els Ship) 


% The acceleration of translation from membrane force in 
the global X coordinate is: 
aximem = 0; 


% The acceleration of translation from membrane force in 
Ene global Yieoeerdinake 1s: 


aYimem = 0; 
% Calculate the angular acceleration about the ship c.g.: 
$ the current contact Dolls: 
CP = [X2 - ((LBP/2)*cos(omega2)), ((LBP/2)- 


L)*sin(omegal)]; 
* So the arm that the force acts through is: 
aa = Sqre (AICA O ae 

* so the angular acceleration due to membrane force is: 
omega dotdotimem - 0; 


* Calculate accelerations due to the Minorsky interaction; 
% The Minorsky force is assumed to act in the direction 


opposite of relative motion. 
% Since this force is in a Gvrrerene, direceion weumucte 


recalculate the virtual mass 


% First, calculate the angle of resultant force 
compared to the ship principal axis 

zetal - wes megal; * the angle of 
Minorsky Lores Eo rhe cStrüucss E EE. 
% Now, calculate the added mass matrix based on this 
angle 

Al = ml* ((all*sim zZz Cal 2 DASS (er 2 
((al1-a22) *sinízetal)*cos (zetal)); ((al1- 
a22) *sin(zetal)*cčos(zetťtal)), (alt*sim( zeta, 2. 
nouos» 
% Combine with physical mass to get the virtual mass 
matrix 


VM1 = M1 + Al; 
* Now calculate the accelerations due to the Minorsky force: 
aXlimin = Fmin*cos (zeta) /VM1 (1,1); 
aYlmin = Fmin*sin(zeta)/VM1(2,2); 
omegal dotdotmin = (82/2/21) "Buin? s io Leto | ee 
omegal) /J66V1; 


% Sum the accelerations due to membrane and Minorsky for the 


portal acelerarionfale 

% due to relative motion and interaction orense men Ec 
axı = aXlmem + aXlmin; 
aYi = aYimem + aYlmin; 


109 


omegal dotdot = omega_dotdotimem + cmegal dotdormin; 


nS DE 
% The acceleration of translation que to membrane roree 
in the globalik eoor aee LS: 
ax2mem = Q0; 
% The acceleration of translation due te menbranee merce 


in the global Y coordinate is: 
aY2umerq - P]:; 


% The angular acceleration about the ship c.g. due to 


membrane force is: 
omega2 dotdotmem = 0; 


% Now calculate accelerations due to the Minorsky 


interaction: 
* The Minorsky force is assumed to act in the direction 


opposite emesclamucm EM 
* Since this force is in a different direction e EN 


recalcula te mier ul mass 


% First, calculate the angle of resultant force 
compared to the ship principal axis 

zeta2 = pi - omega2 + zeta; $ the angle of 
Minorsxy Foco opuestas en gs Da f E 
% Now, calculate the added mass matrix based on this 
angle 

A2 = m2*[(all*sin(zeta2) 2 + az2*ceoc(zZetal eee 
((all-a22)*sin(zeta2)*cos(zeta2)); ((all- 
a22)*siní(zeta2)*cosí(zeta2)), (ali*sin(zeta2) 2 + 
a22*cos(zeta2) 129] 
% Combine with physical mass to get the virtual mass 
mal cios 


VM2 = M2 + A2; 
Y Now calculate the accelerations due to the Minorsky force: 
aX2min = Fmin*cos(zeta-pi) /VM2 (1,1); 
aY2min = Fmin*sin(zeta-pi) /VM2 (2,2) ; 
omega2 dotdotmin = (LBP2/2) *Fmin*sin (omega2- 
Zeta) es 


+ Sum the accelerations due to membrane and Minorsky for the 
total iacolerat on dus 
* due to relative motion and interaction for this time steps 
aX2 = axX2mem « aX2min; 
aY2 = aY2mem + aY2min; 
omega2 dotdot = omega2_ dotdotmem + omega2_dotdotmin; 


% caleta dew ue lo cali 
THI) = E + aXl*step; 
T1(2) = + arlsstep, 


) 
) + aX2*step; 
) + aY2*step; 
omega dotl + omegal dotdot*step; 
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T2302 
DES T 
ee qo 


i 
D 
"aea - 


end 


omega dot2 - omega dot2 + omega2 dotdot*step; 


el ci aereos) =+ reltrans (2) 2)/step; 
Eabs = Eabs + Emin + Emem - Emem last; 
Emem last = Emem; 


Peon aie o See ore cn 
singlecargo 


I 


% single3.m 


5$ Script to perform time-domain analysis for single-hulled 
tanker collis Teon: 

$ The "3" indicates thatithis script iS COrphasce Se S 
con i egon gw ei 

% from the time of impact on the inner longitudinal 
bulkhead, or cargo boundary muta 


* membrane ruptures. 


% Input: all the dynamic variables from srngle2 
% dt from energy m 
% Vi from energy.m 
% V2 from energy.m 
% VM1 from energy.m 
% VM2 from energy.m 
% alpha from vargen.m 


o 


Output: Generated variates and values 


Date created: 3/10/98 
Last updated: 4/23/98 


o? o? 


Determine nearest transverse structures on the inner 


shell, and distance to 
* each for use in applying Van Mater's extension to Jones 


o 


meenod m.: 


o? 


Es Xe 

while BH(j) « L 
¡A 
end 

a = BAHI e I; 

Ir 3 5 A0 
poem T 
else 

[NEM 


end 
* Ensure variable "a" represents the "short leg" of the 


strained plate 


ie s 
cC=ā; 
a=b; 
DEE 

end 


% Calculate deflection at which plate fails, per Van Mater's 


extension to Jones 
defl lim- -0.452*a; 
* Reset deflection counter to zero for thmrsemnevUeNE E SEE 


m 


der — /6- 


% Set new plate thickness for interior bulkhead: 
b uod c MUN OS 


Ao 


6%%% Begin time-step routine %%%%%S%%% 
e abs(defl) « abs(defl lim) 
relvel < endvel 
break 
end 
% Calculate new postions and rotations at end of time step 
time = time + step; 
KI - XT + T(t) *steu; 
<= Yi Steps 
MZ X2 gx q20T)! *step 
po = Y2 = 202) sep 
omegal = omegal + omega dotl*step; 
omega2 = omega2 + omega dot2*step; 
ie = qu Ee 
Pen = Pen + dPen; 
if Pen » Pmax 
Pmax = Pen; 
end 
TE -CLER 
L=LBE> 
break 
end 
ilte 
180. 
break 
end 
if Pen » Beam 
Pen - Beam; 
break 
end 
if Pen « O0 
break 
end 


= Calculate relative translation, perecer c MEM cM a 
impact point in Chis tame seep 


S1 = [T1(1)+((LBP/2) -L) *omega_dotl*sin(omegal), 
T1(2)-((LBP/2) -L) *omega doti*cos(omegal)]; 

S2 = [T2(1) - (LBP2/2)*omega dot2*sin(omega2), 
T2(2) * (LBP2/2) *omega dot2*cos (omega2)] ; 

reltrens, = (92-5) bep 


Ss Calculate direction tor se batt Joe eens owen 
Zeta = aban wees er omeot eel bree 


* Calculate penetration and change in location 
dPen - sqrt((reltrans(1))^2 + 
(reltrans(2))72)*cos(omegal + (3*pi/2) - zeta); 


115 


dL = sqrt((reltrans(1))°2 4+  (reltrans(2)) 2 eco 
+ (3*pi/2) - zeta); 


= Calculate the membraoneweet Veet on 
defl = defl + reltrans(2)+((LBP/2) -L) *sin(omegal) ; 


$ Calculate the resistance forcece O PEN Ene mere 


LE det FEO 
Emem - sigma y*t plate*breadth*defl^2/spacing; 
Fmem = Emem/abs (defl); 
else 
break 
end 
if defl < defl lim 
defl = defl lim; 
Emem - sigma y*t plate*breadth*defl^2/spacing; 
Fmem = Emem/abs (defl); 
rUpcurce2T E 
end 


% Calculate the force resulting from the "Minorsky 
mechanism" 

% Rtt is the total "resistance factor" 

$ dRt is the differential resistance factor for this 
time step 

* depth is the distance of penetration 

* ddepth is the differential distance of penetration 

* t is the aggregate in-plane structure thickness 

% alpha is the bow half-entrance angle 

5 This formula accounts for tne triangular bow wedge 
Cie emer 
and dynamic collision angle, but plic cesna 
limits on striking ship beam. Should be modified so 
that if width exceeds beam, 
remaining area is rectangular.... 


o? oV? 


ov? 


ddepth = sqrt (reltrans(1)^2 - reltrans(2)^2); 

Rtt -(depth^2) *t*tan(alpha*pi/180)/ (1- 
((tan(alpha*pi/180))^2/((tan((omegal-omega2)*pi/180)^2)))); 

depth = depth + ddepth; 

dRt =  (depth^2)*t*tan(alpha*pi/180)/(1- 
((tan(alpha*pi/180))*2/((tan( (omegal-omega2) *pi/180)~2))))+ 
abs (dL) *Pen - Rtt; 


% Calculate the delta-KE from the Minorsky relation; 
Pt) = cee 10 SEM 


t$ The corresponding force is 
Fmin = Emin/abs (ddepth) ; 


Calculate resulting accelerations from the membrane force 
For Phase I, the membrane force is assumed to act 
erpendicularly to the hull surface 
% of the struck ship, and the Minorsky force acts to oppose 
thesdarection of relaci ePmotion: 
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o? oV 


EC) 


tmp e buck ship) 
Calculate the virtual mass 
First, calculate the angle of resultant force 


compared to the ship principal axis 


o? o? ol? 


Zee ou 2. % the membrane force is 
always normal to the struck ship..... 
% Now, calculate the added mass matrix based on this 


angle 

A] - mli*[(all*sin(zetal)^2 4 a22*cos(zetal) 2), 
((a11-a22)*sin(zetal)*cos(zetai)); ((aii- 
a2d2)*etmizetaly *eoo (zeta) ye tall team Zetal ie 2 7 
Ao) * ee Zee 2]. 
E Combine with physical mass to get the virtual mass 


matrix 
VM1 = M1 + Al; 


% The acceleration of translation in the global X 


coordinate is: 
aXlmem = Fmem*sin(omegal)/VM1 (1,1); 


% The acceleration of translation in the global Y 
ese rdinarte mies 
aYlmem = -Fmem*cos (omegal) /VM1 (2,2); 
% Calculate the angular acceleration about the ship c.g.: 
% the current contact point is: 
CP = [X2 - ((LBP/2)*cos(omega2)), ((LBP/2)- 


resin enegal) |i 
* So the arm that the force acts through is: 
arm = sgr- ((XI1-CP(D)) 2 + c YID DOSE 
so the angular acceleration is: 
omega dotdotimem - -((0.5-Lnd)/abs(0.5- 
Lnd))*Fmem*arm/J66V1; 


o? 


* Now calculate accelerations due to the Minorsky 


mac eract on ; 
% The Minorsky force is assumed to act in the direction 


opposite of relative motion. 
% Since this force is in a differ n wee 4 aus: 


recalculate the virtual mass 


% First, calculate the angle of resultant force 
compared to the ship principal axis 

zetal = zeta - omegal; * the angle of 
Minorsky force to the. Sul Do 
% Now, calculate the added mass matrix based on this 
angle 

Al =m*[(alisein(ietall2 + a? costarle 
((al1-a22)*sin(zetal)*cos(zetal)); ((al1- 
a22)*sin(zetal)*cos(zetal)), (ali*sin(zetal) 2 + 
a22*cos(zetal)^2)]; 
% Combine with physical mass to get the virtual mass 
matrix 


VML = ML + AL; 
PIS 


% Now calculate the accelerations due to the Minorsky force: 
aX1min = Fmin*cos (zeta) /VM1(1,1); 
aYlmin = Fmin*sin(zeta) /VM1 (2,2); 
omegal dotdotmin - -((LBP/2)-L)*Fmin*sin(zeta - pi - 
omegal)/J66V1; 


% Sum the accelerations due to membrane and Minorsky for the 


total accleration due 
$ due to relative motion and interaction for this time step 


axl = aXimem + aXlmin; 
aYl = aYlmem + aYlmin; 
omegal dotdot - omega dotdotimem + omegal_dotdotmin; 


Shae 
Calculate the virtual mass for cheo nor sinew ome 
acceleration: 


o? o? 


% First, calculate the angle of resultant force 
compared to the ship principal axis 

zeta2 = omega2 - (omegal + pi/2); 
% Now, calculate the added mass matrix based on this 
angle 

A2 = m2*[(all*sin(zeta2)°2 + a22*ecos(zeta2) 232 
((all-a22)*si1n(zeta2)*cosí(zeta2)) s; (tail - 
a22)*s1in(zeta2) *cos(tetoa2 ea 


al2*cosiuer=a2) 2 


% Combine with physical mass to get the virtual mass 


Moe ra 
VM2 = M2 + A2; 


$ The acceleration of translation in the global X 
coordinate is: 

aX2mem = -Fmem*sin(omegal) /VM2 (1,1) ; 
% The acceleration of translation in the global Y 


coordinate is: 
aY2mem = Fmem*cos (omegal) /VM2 (2,2) ; 


% The angular acceleration about the ship c.g. is: 
omega2 dotdoctmem - (-Piiem*simienegal-eneca = 
pi/2)*(LBP2/2))/J66V2; 


% Now calculate accelerations due to the Minorsky 


interac ilion: 
% The Minorsky force is assumed to act in the direction 


opposite of relative motion. 
* Since this force is in a difrerent Girecrion entes: 


recalculate the virtual mass 


% First, calculate the angle of resultant force 
compared to the ship principal axis 
zeta2 = pl - omega2 + zeta; * the angle of 


Manorsky force to the Strings 


PEG 


% Now, calculate the added mass matrix based on tnis 
angle 
po =em2* (lali*sin(zetaZz) 2 4 a22*cos(zetas2) 2), 
({(ald-a22) *sin(zetaZ) *cos(zeta2)); (lall- 
a: sun (zeta2)*cos (lzeta2)), (all*sin(zetaZ) 2 + 
a22*@os (zetaz) 2)]; 
% Combine with physical mass to get the virtual mass 
maCrix 
VM2 = M2 + A2; 
% Now calculate the accelerations due to the Minorsky force: 
aX2min = Fmin*cos (zeta-pi) /VM2(1,1); 
aY2min = Fmin*sin(zeta-pi1) /VM2 (2,2) ; 
omega2 dotdotmin - (LBP2/2)*Fmin*sin(omega2- 
zeta)/J66V2; 


% Sum the accelerations due to membrane and Minorsky for the 


total accleration due 
* due to relative motion and interaction for this time step 


aX2 = aX2mem + aX2min; 

aY2 = aY2mem + aY2min; 

omega2 dotdot = omega2 dotdotmem + omega2 dotdotmin; 
% Calculate new velocities 

TIU xs TIGE qx ep: 

T3002.) mcd rae YO TET 

VO E 2 qo npo EE. 

T2 (2) = ovg F aY tep, 

omega dot1 - omega dotl + omegal dotdot*step; 


omega dot2 - omega dot2 + omega2 dotdot*step; 
relvel - sqrt(reltrans(1)^2 - reltrans(2)^2)/step; 
Eabs = Eabs + Emin -« Emem - Emem last; 

Emem last - Emem; 


% Determine which tanks have been breached 


Singlecargo 


end 


DI 


% single4.m 


% Script to perform time-domain analysis for single-hulled 


banker col cM 
* The "4" indicates that this script is for phase 4 of the 


Collision welche ts 
% from the time the inner cargo bulkhead ruptures until the 


collision ends. 


* Input: all the dynamic variables from single3.m 
% dt from write.m 

% Va from write.m 

% V2 from write.m 

% VM1 from 

$ VM2 from 

% alpha Peon wiles. ml 


oe 


Output: Generated variates and values 


Date created: 3/10/98 
Last updated: 4/23/98 


oV oe 


$%%%%%S%% Begin time-step routine %%%%%%%%% 
while relvel « endvel 


Fabs - Eabs + Emin; 
* Calculate new postions and rotations at end of time step 
time = time + step; 


Gl = XIT Sie or 


Yio t SURE cece 
X2 = X2 E TUI p 
Vo. = YQ 2-2 (2) Sige 
emegal = omegal + omega doti tete 
omega2 = omega2 + omega _dot2*step; 
D — ea 
Pen = Pen + dPen; 
if Pen > Pmax 
Pmax = Pen; 
end 
CE. cu 
e 
break 
end 
if L » LBP 
LLBP; 
break 
end 


if Pen > Beam 
Pen = Beam; 
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break 
end 
if Pen « O 
break 
end 


* Calculate relative translation, penetration, and change in 
impact point n his time Step 
S Sl 1S He total velocity (from linear and rotational 


modien of Che impact point on Ship i. 
% 9S2 is the same velocity ior Ship 2. 


Si TI s (LBP/2) L)tonega dotigsinlomnsga Nk 
T1 (2) +((LBP/2) -L) *omega_dot1*cos (omegal)] ; 


S2 = [T2(1) - (LBP2/2)*omega_dot2*sin(omega2), 
T2 (2) +(LBP2/2) *omega_dot2*cos (omega2) ] ; 
reltrans = (S2-S1)*step; 


Se Cabkeulateralirecelon Gn relat] Vermrans laracon 
zeta - atan2(reltrans(2),reltrans(1)); 


* Calculate penetration and change in location 
aPen = sqrt ((releraems (i) 2 = 


(reltrans(2))^2)*cos(omegal + (3*pi/2) - zeta) ; 
dL = sgrt((reltrans(1))°2 + (reltrans(2))~ 2) *sin(omegal 
eo ois) 2) ee 


* Calculate the force resulting from the "Minorsky 
mechanism" 
Rtt is the total "resistance factor" 
dRt is the differential resistance factor for this 
time step 

depth is the distance of penetration 

ddepth is the differential distance of penetration 

t is the aggregate in-plane structure thickness 

alpha is the bow half-entrance angle 

This formula accounts for the triangular bow wedge 
Ty; 

and dynamic collision angle, but places no 

limits on striking ship beam: In these cone 
be modified so that if width exceeds beam, 

5$ remaining area is rectangular.... 


o? o? 


geome 


o? o9 (f o9 o? o? o oe 


ddepth = sqrt (reltrans(1)^2 + reltrans(2)2); 

REE = (depth 2) 4b “bam d No ESO AME 
((tan (alpha*pi/180))^2/((tan((omega1-omega2)*pi/180)^2)))); 

depth = depth 4 ddepth; 

GRt = (depth*2) *t*tan(alpha*pi/180) /(1- 
((tan(alpha*pi/180))^2/((tan((omegai-omega2)*pi/180)^2))))- 
abs(dL)*Pen - Rtt; 

* Calculate the delta-KE from the Minorsky relation; 


o 


Ema ardre IO CdR; 
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% The corresponding force is 
Fmin = Emin/abs (ddepth); 


* Ship 1 (Struck ship} 


% The acceleration of translation from membrane force in 
the global X coordinate is: 

axlmem = O5 
% The acceleration of translation from membrane force sme 
the global Y coordinate is: 

aYlmem = 0% 
% Calculate the angular acceleration about the ship Gage 

+ the G@Wrrent Cem@eeee soem 1. . 
CP = [X2 - ((LBP/2)*cos(omega2)), ((LBP/2) - 


L)*sin(omegal)]; 
$ so the arm that the force acts through is: 
arm - SqreL((XI-CPUOIODEZ X1 00 29 7 BN 
5 So the angular acceleration due to membrane force is: 
omega dotdotimem - 0; 


* Calculate accelerations due to the Minorsky interaction; 
% The Minorsky force is assumed to act in the direction 


opposite of relative motion. 
% Since this force 1S in a different direc cio Eye 


recalculate the virtual mass 


% First, calculate the angle of resultant force 
compared to the ship principal axis 

zetal = zeta - omegal; 5$ the angle of 
Minorsky Force to Che sinless D T 
% Now, calculate the added mass matrix based on this 
angle 

Al = ml* [(all*sin(zetal) 2 4 a22*cosi(zetat ae 
((a311-a22)*sin(zetal)*cos(zetal)); “(ail 
a22) *sin(zetal)*cos(zetal1)), (all*sin(zetal)"2 + 
a22 tos (metal 2) |= 
% Combine with physical mass to get the virtual mass 
Mee Cis: 


VMI = ML F Al; 
% Now calculate the accelerations due to the Minorsky force: 


aXimin = Fmin*cos (zeta) /VM1(1,1); 
aYlmin = Fmin*sin(zeta) /VM1 (2,2); 
omegal dotdotmin = -((LBP/2)-L)*Fmin*sin(zeta pE 


omegal)/J66V1; 


% Sum the accelerations due to membrane and Minorsky for the 


total accleration due 
$ due to relative motion and interaction for this time Seen 


axi = aXimem + aXimin; 

aYl = aYimem + aYimin; 

omegal dotdot = omega_dotdotlmem + omegal_dotdotmin; 
% Ship 2 
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% Thoe e llera on OE translation due to membrane force 


in che global xX ‘@@ordinate. 1s: 
aX2mem = O0; 
e —aCee loratton Or translation- due Co membrane force in 


Ene Global Y coordinate is: 

aY2mem = 0; 
% The angular acceleration about the ship c.g. due to 
membrane force is: 

omega2 dotdotmem = 0; 
% Now calculate accelerations due to the Minorsky 
interaction; 
% The Minorsky force is assumed to act in the direction 
opposite of relative motion. 
% Since this force is in a different direction we must 


o 


recalculate the virtual mass 


% First, calculate the angle of resultant force 
compared to the ship principal axis 

zeta2 = pi - omega2 + zeta; 6 the angle of 
Minorsky force EG Ene SEDET E 2. 
% Now, calculate the added mass matrix based on this 
angle 

AQ =] m2* (all*sin(zeta2' 2 + a22 "cos (eran 
((all-a22)*sin(zeta2) *cos(zeta2)); ((al1i- 
a22)*sin(zeta2)*cos(zeta2)), (all*sin(zeta2)"2 + 
a22*cos(zeta2)^2)]; 
% Combine with physical mass to get the virtual mass 
matrix 


MZA AM ee 

% Now calculate the accelerations due to the Minorsky force: 

aX2min = Fmin*cos(zeta-pi)/VM2(1,1); 

av omnes cesta acu UM 

omega2 dotdotmin s (LBP2/2)*Fmin*sin(omega2- 
zeta)/J66V2; 
* Sum the accelerations due to membrane and Minorsky for the 
total accleration due 
* due to relative motion and interaction for this time step 


aX2 = aX2mem + aX2min; 

aY2 = aY2mem + aY2min; 

omega2 dotdot = omega2 dotdotmem + omega2_ dotdotmin; 
% Calculate new velocities 

JI Lue + aXl*step; 

T2 E + aYl*step; 


) 
) + aX2*step; 
27 ) + aY2*step; 

ND Omega coll, - omegal cereer Ep. 
omega dot2'- omega dot? frr omega docdori Een 


1(2 
T2(1) = T2(1 
2(2 


relvel = sqrt(reltrans(1)^2 + reltrans(2)°2)/step; 
Eabs - Eabs + Emin + Emem - Emem last; 
Emem last = Emem; 

% Determine which tanks have been breached 


Singlecargo 


m] 





% Singlecargo.m 


breached during 
% each time step in the simulation for the single hull 


o 


tanker 


% the script to determine which cargo bulkheads have been 


Dace m re O SS 


o? o? 


Last revision: 3/17/98 
zi noue L and Pen from the collision phase scripts 
Quc out: 


o? o 


Flags corresponding to the cargo compartments that will 
release oil 


= 15 
while L > SBH(j) 

Jo Seat; 
end 6 at this end, BH(j) is the first bulkhead aft of 
the damage location 


Jp e» 
Jp, E 
Tank le = 1; 
lf Pen »-24 


du 
end 
end 
end 
IMS — 
3 mupEbEUu5e2 9 
T's li 
end 
end 
if j == 4; 
au EMEND 
Tan. so ME 
end 
end 
Tr] mec 
PEU USO 
Tan qu 
end 
end 
JU ME 
LES RUDLUre =O 
Tank5Sp = I; 
IF rüptüre? m E) 
Tanko T: 
end 
end 
end 


end 


2s i: 
ir erupe rom Ro 
Slop aa 
TE rUpEeUro2 “Seg 
anb le 
end 
end 


% doublei.m 


+= Script to periorm time domain analysis tor double-halzed 


tanker Geol laeaen. 
The "1" indicates that this script is for phase 1 of the 


colle ion Which -is 
% from the time of impact until the shell membrane ruptures. 


o 


o? 


ES. CUE from energy.m 
5 NE Crom energy.m 
% V2 from energy.m 
% VM1 from 
% VM2 from 
6 alpha from vargen.m 


% Output: Generated variates and values 


= Datervereatea- 12/7207 47 
% Last updated: 4/23/98 


* Reset flags for tank breaching to zero: 
TankıP = 
Tankis = 
Tank2P = 
Tank2C = 
Tank3P = 
Tank: ci: 
Tank4C = 
Tank5P.= 
Panes 
FankSP = 
Tank6P = 
mamie = 
ll cM 


". 0589 Ne Ne "mo. We "+ 


we "o 


o mo "a 


OS OO 'DODIDDO DIO o 


o 


% Determine nearest transverse frame structures, and 


distance to each for use in applying 
$ Van Mater’s extension Lo denes met oca 


Feon 

while BH(j) « L 

De eb: 

end 
a BO 
¡O AD 
* Ensure variable "a" represents the "short leg" of the 
strained plate 
SIME a»b 

Cca, 


end 


% Calculate deflection at which plate fails, per Van Mater's 


extension to Jones 
defl lim= -0.452*a; 


% Initialize new variables (Subscript 1 = struck ship, 
subscript 20 Seriki Ship) 

time = 0; 

AT 

Yran 


X2= (LBP/2) -L+(LBP2/2) *cos (phi*pi/180) ; 
Y2= (Beam/2)+(LBP2/2) *sin(phi*pi/180); 
omegal=0; 

omegalao cl = 0; 
omega2=(phi+180)*(pi/180); 

omega dot2 = 0; 


defl = 0; 
= ie ee 
T2 2 Vel2; 


rupture ESOS 
rappeur e E 
rupture3=0; 
Ge 

ddepth = 0; 
dept nsss, 
EOT AER 
Emin = 0; 
Eabs = 0; 
Emem last - 0; 
Ge 

Pen = 0; 
dPen = 0 
Pmax = 0; 


o 


Begin time-step routine %%%%%%%3% 
while abs(defl) < abs{defl lim) 

* Calculate new postions and rotations at end of time step 
time = time + step; 

A SO 

Y1 = Y1 O 

X2 = XZ ee Seco, 

Y2 = Y2 paso seen. 


omegal = omegail. + oméga acti step, 
omega2 = omega2 + omega dot2*step; 
LA yee ca 
Pen = Pen + dPen; 
if Pen > Pmax 
Pmax = Pen; 
end 
if L > LBP 
L-ECBE; 
break 
end 


is Ll. «e 0 
Ie 
break 
end 
if Pen » Beam 
Pen - Beam; 
break 
end 
if Pen « O 
break 
end 


* Calculate relative translation, penetration, and change in 
impact point in this time step 

S1 - [T1(1)«-((LBP/2)-L)*omega doti*sin(omegal), 
T1(2)+((LBP/2) -L) tomega dotl*cos(omegal)]; 


S2 - [T2(1) - (LBP2/2)*omega dot2*sin(omega2), 
T2(2)- (LBP2/2) *omega dot2*cos (omega2)] ; 
reltrans 2 (52991) “step, 


* Calculate direction of relative translation 
zeta = atan2 (reltrans (2) ,reltrans(1)); 


% Calculate penetration and change in location 
dPen = sqrt((reltrans(1))"2 + 


(reltrans(2))^2)*cos(omegal + (3*p1/2) - zeta); 
al. ="sqrt ((reltrans(1)) "2 + (reltrans (2) ) 92) *simteneca 
uec ou RNC SEU 


% Calculate the membrane deflection 
defl = defl + reltrans(2)+((LBP/2) -L) *sin(omegal) ; 


$ Calculate the resistance force of the membrane 

if defl « O 
Emem - sigma y*t plate*breadth*defl 2/spacing; 
Fmem = Emem/abs (defl); 

else 
break 

end 

rswdefl1-«c- desdeslsm 
dele m 
Emem - sigma y*t plate*breadth*defl 2/spacing; 
Fmem - Emem/abs(defl); 
fupture = 21, 

end 


% Calculate the force resulting from the "Minorsky 
mechanism" 

Rtt is the total "resistance factor" 

dRt is the differential resistance factor for this 
time step 

depth is the distance of penetration 

ddepth is the differential distance of penetration 
t is the aggregate in-plane structure thickness 
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o? o? 


o9 o9 o? 


pha is the bow half-entrance angle 
his formula accounts tor theta Ulead e 


nd dynamic collision angle, but places no 
imits on striking ship beam. Should be modified so 
that 1f width exceeds bean, 

% remaining area is rectangular... 1 


al 
T 
geometry 
an 
E 


ddepth = sqrt (reltrans(1)"2 + reltrans(2)*2); 

Rtt = depo ud ee 
((tCan(alpha*pi/180))^2/((tan((omegal-omega2)*pi/180)^2)))); 

depth = depth + ddepth; 

dRt = Ce risu uan apa 
((tan(alpha*pi/180))^2/((tan((omegal-omega2)*pi/180)^2))))- 
abs (dL) *Pen - Rtt; 


% Calculate the delta-KE from the Minorsky relation; 
Emin —- dNE*IO 6*dEEC 


% The corresponding force is 
Fmin = Emin/abs(adepth); 


Calculate resulting accelerations from the membrane force 
For Phase I, the membrane force 1S assumed to act 
erpendicularly to the hull surface 
of the struck ship, and the Minorsky force acts PONES 
the direction of relative motion. 


o? o? 


oe ! (Cj 


np EE Ium) 
Calculate the virtual mass 
First, calculate the angle of resultant force 
compared to cne ship Piane palmals 


o? o? o? 


Zeta ls = way 2 % the membrane force is 
always normal to the struck ship..... 
% Now, calculate the added mass matrix based on this 
angle 

Al - mi*[(alli*sin(zetal)^2 + a@2*cos(Zetal), Ze 
((all-a22) *sin(zetal)*cos(zetal)); ((all- 
a22)*sin(zetal)*cos(zetal)), (ali*sin(zetal)"2 + 
az2* Cos (zetal) 2) i 
% Combine with physical mass to get the virtual mass 
Wena as 

VM1 Ml + Al; 
% The acceleration of translation in the global X 


coordinate is: 
aXlmem = Fmem*sin(omegal) /VM1 (1,1); 


$ The acceleration of translation in the global Y 
coordinate is: 

aYlmem = -Fmem*cos (omegal) /VM1 (2,2) ; 
% Calculate the angular acceleration about the ship Cc E 


*v-the current 'conbuet sed uS 
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CP - [X2 - ((LBP/2)*cos(omega2)), ((LBP/2)- 
L)*sin(omegal)]; 
* so the arm that the force acts through is: 
arm = scc (X1 CP(1)) 2-+ (yq - Gp(29892); 
% so the angular acceleration is: 
omega dotdotimem - -((0.5-Lnd)/abs(0.5- 
Lnd))*Fmem*arm/J66V1; 


% Now calculate accelerations due to the Minorsky 


DIEI iro; 
* The Minorsky force is assumed to act in the direction 


opposite of relative motion. 
< Since this force ds $2nm'wdlitferentMirectio MES 


recalculate the virtual mass 


% First, calculate the angle of resultant force 
compared to the ship principal axis 

zetal = zeta - omegal; * the angle of 
morsky roce eo tner” steel ships 
% Now, calculate the added mass matrix based on this 
angle 

A1 - mi*[(all*sin(zetal) 2 + a22*cos(zetal) 2), 
C(all-a22)*sin(zetal)*cos(zetau9))- ut agu- 
a22)*sin(zetal)*cos (zetal)), (all cin Z-ca Je 
a22*cos(zetal) 2) |; 
% Combine with physical mass to get the virtual mass 
matrix 


IIS a 
* Now calculate the accelerations due to the Minorsky force: 
aXlmin = Fmin*cos (zeta) /VM1(1,1) ; 
aYlmin = Fmin*sin(zeta) /VM1 (2,2); 
omegal dotdotmin = -((LBP/2)-L)*Fmin*sin(zeta - pi - 
omegal)/J66V1; 


* Sum the accelerations due to membrane and Minorsky for the 


total accleration due 
$ due to relative motion anda interaction for thisWEIme Ste 


axl = aXimem « aXimin; 

aYl = aYlmem + aYimin; 

omegal dotdot - omega dotdotimem « omegal dotdotmin; 
$5 9 59 6 0 6 6 6 66 6 66 6 5 6 6 56 6 0 66 6 66 6 5 6 6 5 6 6 0 6 6 0 0 6 6 5 6 6 6 0 
C TTE 
% 


Calculate the virtual mass for the membrane force 
acceleration: 


% First, calculate the angle of resultant force 
compared to the ship principal axis 

zeta2 = omega2 - (omegal + pi/2); 
$ Now, calculate the added mass matrix based on this 
angle 

A2 = m2*{(all*sin(zeta2) "2 + a22*cos(zeta2) 2), 
((al1l-a22) *sin(zeta2) *cos(zeta2)); ((all- 
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a22)*sin'(zeta2)*cos(zeta2) (aie ea eee ee 
awe * emi seas 


% Combine with physical mass to get the virtual mass 


matrix 
VM2 = M2 + A2; 


% The acceleration of translation in theme lesa. 
coordinate is: 

aX2mem = -Fmem*sin(omegal) /VM2 (1,1); 
% The acceleration of translation in the global Y 


coordinate is: 
aY2mem = Fmem*cos (omegal) /VM2 (2, 2) ; 


The angular acceleration about the ship c.g. is: 
omega2 dotdotmem = (-Fmem*sin(omegal-omega2- 
¡A EP 27 2 EN occ ve. 


no 


% Now calculate accelerations due to the Minorsky 


interaction, 
% The Minorsky force is assumed to act in the direction 


opposite of relative motion. 
% Since this force 1s 1n a ditiérent direction. ence 


recalculate the virtual mass 


First, calculate the angle of resultant force 
compared to the ship principal axis 

zeta2 = pl - omega2 + zeta; * the angle of 
Minorsky force to the striking o hip eee 
% Now, calculate the added mass matrix based on this 
angle 

A2 = m2*[(all*sin(zeta2) 2 + a22*cos (zeta se 
(dall a22): Sin Ze ea aoe l zea a EET a 
a22)*sin(zeta2) *cos(zeta2)), (all*sin(zeta2) 2 + 
Se2272GO> (Zonas, Jia 
% Combine with physical mass to get the virtual mass 
matr rk 


MNA MSIE 
* Now calculate the accelerations due to the Minorsky force: 
ax2mrm = BEntiu*tcostzesd-91)/VM2WN 
aY2min = Fmin*sin(zeta-pi) /VM2 (2,2) ; 
omega2 dotdotmin = (LBP2/2) *Fmin*sin (omega2- 
zeta) /J66V2; 


% Sum the accelerations due to membrane and Minorsky for the 


total accleration due 
* due to relative motion and interaction for this time Steg 


aX2 = aX2mem + axX2min; 

aY2 = aY2mem + aY2min; 

omega2_ dotdot = omega2_ dotdotmem + omega2_dotdotmin; 
$ Calculate new velocities 

TL) “TILA | Step 

TL(2) = T1(2) + aa Seer, 
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ToT) T2 (1) + aX2*step; 

pL = T22). + aY2* step; 

omega dotl - omega dotl1 + omegal_dotdot*step; 
Occ o E-Emnedgaeeoto c omnegazecobdot*step: 
relvel - sqrt(reltrans(1)^2 « reltrans(2)^2)/step; 
Eabs - Eabs + Emin + Emem - Emem last; 

Emem last - Emem; p 


% Determine which tanks have been breached 
doublecargo 


end 


% double2.m 


Script to perform time-domain analysis for double-hulled 


tanker collision. 
The "2" indicates that this script ISmo oree ee 


col liston wile a is 
% from the time the outer shell membrane ruptures until the 


o 


inner bulkhead is conta ertecas 


o 


% Input: all the dynamic variables from doublel.m 
$ dt from writem 

% V1 from write m 

% va from write.m 

% VM1 Trom 

% VM2 Erom 

% alpha from write.m 


% Output: Generated variates and values 


% Date created: 11/3/97 
* Last updated: 4/23/98 


o 
ov 
o? 
o? 
o? 
o 
o\? 
o\? 
o? 


$%%%%S%%% Begin time-step routine 
Pen « DS1 

lvel « endvel 

break 


% Calculate new postions and rotations at end of time step 
time = time + step; 

ie Me Pe) Seo: 

YI Y1 + T1(2)*step; 

XO a) A S 

Y2. =" Yee 922) *sren- 


I 


omegal = omegal + omega dotl*step; 
omega2 = omega2 + omega dot2*step; 
i = TL Sci 
Pen = Pen + dPen; 
if Pen » Pmax 
Pmax = Pen; 
end 
Jub < 20 
EC 
break 
end 
if L » LBP 
LLBP; 


break 

end 

if Pen > Beam 
Pen = Becam, 


break 
end 
if Pen < 0 
break 
end 


% Calculate relative translation, penetration, and change in 
impact point in this time step 

$151.18 the total velocity (Irom linear and Eo cue 
motion) Ob Co- eana Polne On Sato 1s 

* S2 is the same velocity for Ship 2. 


S1 = [T1(1)+((LBP/2) -L) *omega_dotl*sin(omegal), 
T1 (2) +((LBP/2)=L) *omega_deti*cos (omegald'; 


S2 - [T2(1) - (LBP2/2)*omega dot2*sin(omega2), 
T2(2)-(LBP2/2) *omega dot2*cos (omega2)]; 
sel bras =o see aS ele, 


3% Calculate direction of relative translation 
zeta = atan2(reltrans(2),reltrans(1)); 


% Calculate penetration and change in location 
dPen = sqrt((reltrans(1))"2 + 


(reltrans(2))%72)*cos(omegal + (3*pi/2) - zeta); 
dL = sgqrt((reltrans(1)) "2.2 (reletrans (21052) mmu sm ll 
Goi) e iQ 


% Calculate the force resulting from the "Minorsky 
mechanism" 

Rtt is the total "resistance factor" 

dRt is the differential resistance factor for this 


time step 

< depth rs the droan c- Tor Pe ea IN 

* ddepth is the differential distance of penetration 

* t is the aggregate in-plane structure thickness 

* alpha is the bow half-entrance angle 

* This formula accounts for the triangular bow wedge 
geometry 

* and dynamic collision angle, but places no 

* limits on striking ship beam. In chewelutime could 


be modified so that if width exceeds bean, 
6 remaining area 1s rectangular.... 


ddepoth = sqrt (reltrans (1) 2a rele rama 2 eee 

Rtt = (depth^2)*t*tan(alpha*pi/180)/(1i- 
((tan(alpha*pi/180))^2/((tan((omegal-omega2)*pi/180)^2)))); 

depth = depth + ddepth; 

dRt =  (depth^2)*t*tan(alpha*pi/180)/(1- 
((tan(alpha*pi/180) ) *2/((tan( (omegal-omega2) *pi/180)~2))))+ 
abs (dL) *Pen - Rtt; 
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$ Calculate the delta-KE from the Minorsky relation; 
Emin - dKE*10^6*dRt; 


* The corresponding force is 
Fmin = Emin/abs(ddepth) ; 


+ Ship 1 (Struekmshi pe) 


% The acceleration of translation from membrane force in 
the global X coordinate is: 

aximem = 0; 
% The acceleration of translation from membrane force in 
the global Y coordinates. 

aYlmem = 0; 
% Calculate the angular acceleration about the ship c.g.: 

+ the cusmsent Cenbace some S 
CP - [X2 - ((BBP/2)*cos(omega2 BE o 


L)*sin(omegal)]; 

* so the arm that the force acts through is: 
arm = Sqrt ((xXP=CP(1))°2 yt cuum» E 

% so the angular acceleration due to membrane force is: 
omega dotdotimem = 0; 


% Calculate accelerations due to the Minorsky interaction; 
6 The Minorsky force ris assumed to ae In@the directrkam 
opposite of relative motion. 

* Since this force is in a different direction we must 
recalculate the virtual mass 


% First, calculate the angle of resultant force 
compared to the ship principal axis 

zetal = zeta - omegal; “ the angle of 
Minorsky force to the struck ship. wm... 
% Now, calculate the added mass matrix based on this 
angle 

Al = mi*[(all*sin(zetal) “2 7 a22*cos(zeťtal) 20S 
((a11-a22)*sin(zetal)*cos(zetal)); ((aii- 
a22)*siní(zetal)*cos(zetal1)), (ali*sin(zetal) 2 + 
a22*cos(zetal) 2)]; 
% Combine with physical mass to get the virtual mass 
matrix 


VM1 = M1 + Al; 
% Now calculate the accelerations due to the Minorsky force: 
aXimin = Fmin*cos (zeta) /VMi (1,1); 
avYimin = Emin*sin(zeta) WWMi(2 2) 
dim dot dora ia = (LBP/2)-L)*Fmin*sin(zeta - pi - 
omegal) /J66V1; 


* Sum the accelerations due to membrane and Minorsky for the 


total accleration due 
* due to relative motion and interaction for this time Suenan 


axl = axlmem + aximin;, 
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aYl = aYlmem + aYlmin; 


omegal dotdot - omega dotdotimem + omegal dotdotmin; 

CE MEE 
zl The acceleration of translation due to membrane force 
in the global X coordinate is: 

ax2mem = O0; 
% The acceleration of translation due to membrane force 
in the global Y coordinate is: 

aY2mem = O0; 
% The angular acceleration about the ship c.g. due to 
membrane force is: 

omega2 dotdotmem = 0; 


* Now calculate accelerations due to the Minorsky 


interaction. 
% The Minorsky force is assumed to act in the direction 


opposite of relative motion. 
% Since this force is in a different direction we must 


recalculate the virtual mass 


% First, calculate the angle of resultant force 
compared to the ship principal axis 

zeta2 = pl - omega2 + zeta; % the angle of 
bonorskewv rorce to ehe “St mica cto e 
% Now, calculate the added mass matrix based on this 
angle 

A2 =m2* [(ali*sin(zeta2) 2 + ar? cost zeua ee 
((a11-a22)*sin(zeta2)*cos(zeta2)); ((a1i1- 
a22)*sin(zeta2)*cos(zeta2)), (all*sin(zeta2)2 + 
a22*cos(zeta2) 2)1]; 
% Combine with physical mass to get the virtual mass 
Matrix 


VM2 = M2 + A2; 
* Now calculate the accelerations due to the Minorsky force: 
aX2min = Fmin*cos(zeta-pi) /VM2(1,1); 
aY2min = Fmin*sin(zeta-pi) /VM2 (2,2); 
omega2 dotdotmin = (LBP2/2) *Fmin*sin (omega2- 
zeta) /J66V2; 


* Sum the accelerations due to membrane and Minorsky for the 


total accleration due 
* due to relative motion and interaction for this time step 


aX2 = aX2mem + aX2min; 

aY2 = aY2mem + aY2min; 

omega2 dotdot = omega2 dotdotmem + omega2 dotdotmin; 
% Calculate new velocities 

TIT IAE NE EN STE 

El (2) -— quay MS UTOR 

Iq 0g cR LESTIE 

"Vos E T aD ES 


JE» 


omega doti - omega dotl + omeégal detdet*step,; 
omega dot2 - omega dot2 -« omega2 dotdot*step; 


relvel = sqrt (reéltrans(1) 2 4 velerams (2) (7 Sco 
Babs = Fabs + Emin + Emem — Bmem face, 
Emem last = Emem; 


* Determine which tanks have been breached 
doublecargo 


end 


% double3.m 


Script to perform time-domain analysis for single-hulled 


tanker collision. 
% The "3" indicates that this script is for phase 3 of the 


COllacienme When le 
$ from the time of impact on the inner shell, or cargo 


boundary until that 
6 membrane ruptures. 


Input: all the dynamic variables from double2.m 


% dt from energy.m 
% VI from energy.m 
% V2 from energy.m 
% VM1 from energy.m 
% VM2 from energy.m 
% alpha from vargen.m 


Output: Generated variates and values 


o 


Date created: 12/30/97 
Last updated: 4/23/98 


o? o 


% Determine nearest transverse structures on the inner 
shell, and distance to 

5$ each for use in applying Van Mater's extension to Jones 
method.... 


ME 
while BH(j) « L 
J = J+l; 
end 
ae DE S, 
Wr. Ee 
b = L = BH(j-21); 
else 
ee 


end 
¢ Ensure variable "a" represents the "short leg" of the 


strained plate 


loa a>b 
C=a; 
a=b; 
b=c; 
end 


* Calculate deflection at which plate fails, per Van Mater's 


extension to Jones 
defi lim-T 0 doia 


% Reset deflection counter to zero for tals mer mmena=ner 
Get] Ser 


o? 


$%6%% Begin time-step routine %%%%%%%% 
while abs(defl) < abs(defl lim) 
f relvel « endvel 
break 
end 
% Calculate new postions and rotations at end of time step 
time = time + step; 
X1 - X1 + T1(1)*step; 
Y= Ye ri ED 
K2 X2 cp DOS TEE or 
Y2 = Y2 + T EE 
omegal - omegal « omega doti*step; 
omega2 - omega2 + omega dot2*step; 
Door INS 
Pen = Pen + cPenc 
if Pen > Pmax 
Pmax = Pen; 
end 
if L > LBP 
LLBP; 
break 
end 
IE ec Aen 
(0) 5 
break 
end 
1f Pen > Beam 
Pen = Beam; 
break 
end 
if Pen < 0 
break 
end 


$ Calculate relative translation, penetration, and change in 
impact point in this time step 

Sl = [T1(1)+((LBP/2) -L) *omega_dot1l*sin(omegal), 
T1(2)+((LBP/2)-L) *omega dot1*cos (omegal) J; 


S2 = [T2(1) - (LBP2/2)*omega dot2*sin(omega2), 
T2(2)+(LBP2/2) *omega_dot2*cos (omega2)]; 
re leas Pol Seer, 


% Calculate direction of relative translation 
zeta - atanZ2(reltrans(2),reltrans(1)); 


$ Calculate penetration and change in location 
dPen = sqrt((reltrans(1))"2 + 


(reltrans(2))°2)*cos(omegal + (3*pi/2) - zeta); 
dL = sqrt ((reltrans(1))°~2 4 (reltrans(2)) 2)*sin(omeq am 
+ (3*pi/2) - zeta); 


= Caletbare the membrane deflecticn 
defl = defl + reltrans(2)+((LBP/2) -L) *sin(omegal) ; 


~ Calenlavewene sees!  seaomce fonce cot rhesmembsare 


2 feti 0 
Emem - sigma y*t plate*breadth*defl^2/spacing; 
Fmem = Emem/abs (defl); 
else 
break 
end 
if defl < defl lim 
defl = defl him, 
Emem - sigma y*t plate*breadth*defl^2/spacing; 
Fmem = Emem/abs (defl); 
lec 
erd 
$ Calculate the force resulting from the "Minorsky 
mechanism" 


REE tS ther coca restceonece =actor” 
dRt sS the difrerential resistance factorii on ie 


cime step 

% depth is the distance of penetration 

% ddepth is the differential distance of penetration 

% t is the aggregate in-plane structure thickness 

% alpha is the bow half-entrance angle 

% This formula accounts for the triangular bow wedge 
geometry 

% and dynamic collision angle, but places no 

% limits on striking ship beam. Should be modified so 
that if width exceeds beam, 


No pu. 


remaining area is rectangular.... 


ddepth = sqrt (réeltrvans(1) 2 + seltrans(2 2 
Rtt - (depth^2)*t*tan(alpha*pi/180)/(1i- 
((tan(alpha*pi/180))72/((tan((phi+omegal- 
emegaZ) *p1/meoyr 2))))% 
depth = depth + ddepth; 


GaRt =  (depth^2)*t*tan(alpha*pi/180)/(1- 
((tan(alpha*pi/180) )*2/((tan((phi+omegai- 
omega2)*pi/180)^2))))-« abs(dL)*Pen - Rtt; 


* Calculate the delta-KE from the Minorsky relation; 
Ent ES CGU 


% The corresponding force is 
emin —- Ein, aoe eM 


Calculate resulting accelerations from the membrane force 
For Phase I, the membrane force 1s assumed to act 
erpendicularly to the hull surface 

of the struck ship, and the Minorsky forcera cts Oo OPeee— 
the direction of relative motion. 


o? oo 


o FJ 


So Se tle ato 


Calculate the virtual mass 
First, calculate the angle of resultant force 


compared to the ship principal axis 


o? oe 


Zebal = oly 2 % the membrane force is 
always normal to the struck ship..... 
% Now, calculate the added mass matrix based on this 
angle 

Al = mi*[(alli*sin(zetal)^2 « a22*cos(zeta1)^2), 
((all-a22) *sin(zetal) *cos(zetal)); ((all- 
a22)*sin(zétal) *ees(zetal)), (all*sini rec M 
a22*cos(zetal)^2)]; 
% Combine with physical mass to get the virtual mass 
Mater ie 

VM1 = M1 + Al; 
% The acceleration of translation in the global x 


coordinate is: 
aXlmem = Fmem*sin(omegal) /VM1 (1,1) ; 


% The acceleration of translation in the global Y 
coordinate is: 
aYlmem = -~Fmem*cos (omegal) /VM1 (2,2); 


Calculate the angular acceleration about the ship cra 
* the current contact point is: 
CP = [X2 -"X(LBPZ2)*cosvtomedga2 DNE EPA E 
L)*sin(omegal)]; 
* So the arm that the force acts through is: 
art = SQre((X1l=CP (bh) ) 27-3 (v1 o Ce) 7 re 
% so the angular acceleration is: 
omega dotdotimem = -((0.5-Lnd)/abs(0.5- 
Ind) )*Fmem*arm/J66V1; 


o? 


$ Now calculate accelerations due to the Minorsky 


interaction; 
The Minorsky force is assumed to act in the direction 


opposite of relative motion. 
* Since this force is in a different direction we muse 


recalculate the virtual mass 


o H- 


% First, calculate the angle of resultant force 
compared to Ene ship: prencipalyacis 

zetal = zeta - omegal; * the angle of 
Minorsky Fores tO Ehe steuck ship. 2a 
% Now, calculate the added mass matrix based on this 
angle 

Al - ml*[(all*sin(zeta1)^2 + a22*cos(zetal) 2), 
((all-a22)*sim(zetal)*cos(zeta3)); 9 S07 
a22)*sin(zetal)*cos(zetal)), (all*sin(zetal) 2 - 
a22*cos(zetai) 2)]; 
% Combine with physical mass to get the virtual mass 
maeri 


VM1 = M1 + Al; 
% Now calculate the accelerations due to the Minorsky forece 
aXimin = Fmin*cos(zeta) /VM1 (1,1); 
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Edd Fmin*Sin(zeta) /VM1(2,2) ; 
onecoa Edo doc man LD =P mee aes 
omegal)/J66V1; 


$ Sum the accelerations due to membrane and Minorsky for the 


total accleration due 
* due to relative motion and interaction for this time step 


aX1 - axXlmem + aXlmin; 

aYl = avYigpem + aYimin; 

omegal dotdot - omega dotdotimem « omegal dotdotmin; 
SD 


o oV? 


Calculate the virtual mass for the membrane force 
acceleration: 


$ First, calculate the angle of resultant force 
compared to the ship principal axis 

zeta2 = omega2 - (omegal + pi/2); 
% Now, calculate the added mass matrix based on this 
angle 

A2 = m2*[(all*sin(zeta2)"2 + a22*cos(zeta2) 2), 
((a11-a22)*sin(zeta2)*cos(zeta2)); ({(all- 
a22)*sin(zeta2)*cos(zeta2)), (all*sin(zeta2) "2 + 


pec wet M 


% Combine with physical mass to get the virtual mass 


matrix 
VM2 = M2 + A2; 


% The acceleration of translation in the global X 
coordinate is: 
aX2mem = -Fmem*sin(omegal)/VM2 (1,1); 
% The acceleration of translation in the global Y 
Geordinate iS: 
aY2mem = Fmem*cos (omegal) /VM2 (2, 2) ; 
% The angular acceleration about the ship c.g. 1s: 
omega2 dotdotmem = (-Fmem*sin(omegal-omega2- 


poa (Bee USC E 


% Now calculate accelerations due to the Minorsky 


interaction; 
o : : ; : : 
* The Minorsky force is assumed to act in the direction 


opposite of relative motion. 
5$ Since this force is in a different direction we must 


recelenlate vette ura mass 


6 First, calculate the angle of resultant force 
compared to the ship principal axis 

zeta2 = pi - omega2 + zeta; $ the angle of 
Minorsky foro bo tuesscmISNIgTS)D aE 
% Now, calculate the added mass matrix based on this 
angle 
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A2  m2*[(all*sin(zeta2) 2 «4 a22*cos(zeta2) 2), 
((al1=3422) *sintzeta2) "Costera DN 
a22)*sin(zeta2)*cos(zeta2)), (ali*sin(zeta2) 2 + 
a22*cos(zeta2)^2)]; 

% Combine with physical mass to get the virtual mass 
lee i 

VM2 = M2 + A2; 

* Now calculate the accelerations due to the Minorsky force: 
aX2min = Fmin*cos (zeta-pi) /VM2(1,1) ; 
ay2min -~ Emisora palo M M ee 
omega2 dotdotmin = (LBP2/2) *Fmin*sin (omega2- 

zeta) /J66V2; 


% Sum the accelerations due to membrane and Minorsky for the 
totalo Cler Io Gie 
% due to relative motion and interaction for this time step 
aX2 = aX2mem + aX2min; 
aY2 = aY2mem + aY2min; 
omega2 dotdot = omega2 dotdotmem + omega2 dotdotmin; 


% Calculate new velocities 
Trop o O o ao 
Ti (2) 92) Fi 2) eee Seer. 
T2418) UE) PARAS ESO 
T2442) ea -T242) 4 aves Stee, 


omego aot omega doti + someda lade: dor setae. 
omega dot2 = omega dot2 + omegaZ dotdoe*step, 


relvel = sqrt (reltrans(1) 2.7 reltrans(z F 
Eabs = Eabs + Emin + Emem - Emem last, 
Emem last = Emem; 


% Determine which tanks have been breached 
doublecargo 


end 


% double4.m 


* Script to perform time-domain analysis for double-hulled 


o 


tanker collision. 
% The "4™ indicates that this script is for phase 4 olithe 


collision whichis 
% from the time the inner shell membrane ruptures until the 


o 


inner cargo bulkhead is 
* contacted. 


% Input: all the dynamic variables from double3.m 
% dt from write.m 

% V1 from write.m 

% N2 from write.m 

% VM1 from 

$ VM2 Erom 

% alpha from write.m 


Output: Generated variates and values 


ov? 


Date created: 1/5/98 
Last updated: 4/23/98 


ae A 


$%%%%%% Begin time-step routine %%3%%%%%%% 
while Pen « DS2 
f relvel « endvel 


a 
break 
end 
Fabs = Eabs + Emin; 


* Calculate new postions and rotations at end of time step 
time = time + step; 
A= Xd sse em 
Yl = Yl + T1(2)*step; 
oco a seer 
y2 =) Y204-42(2) “step, 
omegal = omegal -« omega dotl*step; 
omega2 - omega2 « omega dot2*step; 
E 
Pen = Pen + dPen; 
if Pen » Pmax 
Pmax - Pen; 
end 
LE UE E 
[eos 
break 
end 
if L > LBP 
LLBP; 
beak 


end 
1£ Pen > Beam 


Pen = Beam; 
break 
end 
if Pen « O 
break 
end 


% Calculate relative translation, penetration, and change in 
impact point in this time step 
% S1 is the total velocity (from linear and rotational 


motron) of thesmpascbpolmnsec SH E 
% S2 is the Same velocity term shies. 


S1 = [T1(1)+((LBP/2) -L) *omega_doti*sin(omegal), 
T1(2)-((LBP/2)-L) *omega doti*cos(omegal)]; 


S2 = [T2(1) - (LBP2/2)*omega dot2*sin(omega2), 
T2 (2) +(LBP2/2) *omega_dot2*cos (omega2)]; 
reltrans = (S2-S1)*step; 


* Calculate direction of relative translation 
zeta - atan2(reltrans(2),reltrans(1)); 


% Calculate penetration and change in location 


dPen = sqrt ((rellerane(i)) 245 
(reltrans(2))72)*cos(omegal + (3*pi/2) - zeta); 

dL = sqrt((reltrans(1))°2 + (reltrans(2))~2)*sin(omegam 
+ (3*pi/2) - zeta); 


$ Calculate the force resulting from the "Minorsky 
mechanism" 
Rtt is the total "resistance factor" 
dRt is the differential resistance factor for this 
time step 

depth is the distance of penetration 

ddepth is the differential distance of penetration 

t is the aggregate in-plane structure thickness 

alpha is the bow half-entrance angle 

This formula accounts for the triangular bow wedge 
ry 

and dynamic collision angle, but places no 

limits on striking ship beam. In the future, comia 
be modified so that if width exceeds beam, 

* remaining area is rectangular.... 


o? o 


geome 


o? oV (f o? o? o? o9 o? 


ddepth = sqrt (reltrans(1)°2 + reltrans(2) 2) ; 

Rtt = (depth*2) *t*tan(alpha*pi/1i80) / (1- 
( (tan (alpha*pi/180))7*2/((tan( (omegal-omega2) *pi/180)~2)))); 

depth = depth + ddepth; 

dRt = (depth*2) *t*tan(alpha*pi/180) / (1- 
((tan(alpha*pi/180))*2/((tan( (omegal-omega2) *pi/1i80) ~2))))+ 
abs (dL) *Pen - Rtt; 
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% Calculate the delta-KE from the Minorsky relation; 
Emin - dKE*10^6*dRt; 


% The corresponding force is 
Fmin = Emin/abs (ddepth) ; 


= Shigeemisenuck Ship) 


% The acceleration of translation from membrane force in 
the global X coordinate is: 

aXlmem = 0; 
$ The acceleration of translation from membrane force in 
the global Y coordinate is: 

ayimem = 0; 
$ Calculate the angular acceleration about the ship c.g.: 

vehe Curren cone aci DO Pills ocr 
CP = [X2 <= I(LBP/2)+*c0s (omegaz 1 Mee e 


L)*sin(omegal)]l; 
% so the arm that the force acts through is: 
azm = Sqre ((Xi=CP(ij) 92 2 2 yg = Cea E 
% so the angular acceleration due to membrane force is: 
omega dotdotimem = 0; 


% Calculate accelerations due to the Minorsky interaction; 
% The Minorsky force is assumed to act in the direction 


opposite of relative motion. 
% Since this force as in a different direction wa must 


recalculate the virtual mass 


% First, calculate the angle of resultant force 
compared to the ship principal axis 

zetal = zeta - omegal; % the angle of 
Minorsky torce/to' the seee Si 
% Now, calculate the added mass matrix based on this 
angle 

AL = mi*[((a11*Ssia zeta) 2 22274665 (| Zeea) 2 
((al1-a22) *sin(zetal)*cos[(zetal)); ((al1- 
a22)*sin(zetal)*cos(zetal)), (ali*sin(í(zetal1)"2 + 
d22dmco e EE 
% Combine with physical mass to get the virtual mass 
matrix 


VML- =- MI TALT: 
% Now calculate the accelerations due to the Minorsky force: 
aXlmin = Fmin*cos (zeta) /VM1(1,1); 
avimin = FEmin*sin (zeta VM mF 
omegal dotdotmin - (LBP/2)-L)*Fmin*sin(zeta - pi - 
omegal)/J66V1; 


% Sum the accelerations due to membrane and Minorsky for the 


o 


total acelerar ion due 
% due to relative motion and interaction for this time step 
axl = aXlmem + aXlmin; 
aYl = aYlmem t arimin: 
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omegal dotdot - omega dotdotimem « omegal dotdotmin; 


SS 
% The acceleration of translation due to membrane force 


in the global X coordinate is: 
aX2mem = 0; 


% The acceleration of translation due to membrane force 
in gene global Y coordi mares: 

aY2mem = 0; 
% The angular acceleration abeæeut the ship c-g. dueto 
membrane force is: 

omega2 dotdotmem = 0; 


% Now calculate accelerations due to the Minorsky 


interaction; 
o 8 * , . * 
% The Minorsky force ıs assumed to act in the direction 


opposite of relative motion. 
% Since this force 16 In a dilftferemtadgeccrore Am 


recalculate ci a Ire Ua amie se 


% First, calculate the angle of resultant force 
compared to the ship principal axis 

zeta2 = pi - omega2 + zeta; % the angle mai 
Minorsiv c torcestosthessbur mise sim A 
% Now, calculate the added mass matrix based on this 
angle 

A2 z m2*[(all*sin(zeta2)2 «4 a22*cos(zeta2) 2), 
((al1-a22)*sin(ze-a2)*cos(zeta2)); ((a11- 
a22)*sin(zeta299cos(zeta2)), (lall*sinizetaz e 
a22*žcos (zeta2)^2)]l; 
% Combine with physical mass to get the virtual mass 
matrix 


VM2 STMI EAE 
% Now calculate the accelerations due to the Minorsky force: 
ax2min = Fmin*cos(zeta-pi)/ VM ieee 
aY2min = Fuste po 2E 
omega2 dotdotmin = (LBP2/2) *Fmin*sin (omega2- 
zeta) /J66V2; 


% Sum the accelerations due to membrane and Minorsky for the 
total accleration due 
* due to relative motion and interaction for this time step 
aX2 = ax2mem + aX2min; 
aY2 = aY2mem + aY2min; 
omega2 dotdot = omega2_ dotdotmem + omega2_dotdotmin; 


S Calculate new velocities 
TI(Cb) c TIMO?  ax-Estep. 
Ti (2) "2 99) 4 are 
T2(1) s T2(Q1) - oca» sep 
T2112). m T2 (2) ar rS e pP 


omega doti $$ omegadodotdep Scop 
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omega dot 


omega dot2 - omega dot2 + omega2 dotdot*step; 


DR le sooner ce 057 reltranst2) 2 seb. 
Eabs = Eabs + Emin + Emem - Emem last; 
Emem last - Emem; 


$ Determine which tanks have been breached 
doublecargo 


end 
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% double5.m 


* Script to perform time-domain analysis for single-hulled 


tanker collision. 
* The "5" indicates that this script is for phas eer ee as 


collisuonuwhrehgws 
4 from the time of impact on the inner cargo bulkhead, 


Unas 
* membrane ruptures. 


Input: all the dynamic variables from double4.m 


* dt from energy.m 
% E from energy.m 
% V2 from energy.m 
% VM1 from energy.m 
% VM2 from energy.m 
6 alpha from vargen.m 


o? 


Output: Generated variates and values 


Date created: 3/10/98 
Last updated: 4/23/98 


lid 


* Determine nearest transverse structures on the shell, and 


distance to 
* each for use in applying Van Mater's extension to Jones 


Que uc = 


= L; 
while BH(j) < L 
Je o irl; 
end 
a = BHE- L; 
IE Jo 
b = L - BH(j-1); 
else 
DSO, 


end 
% Ensure variable "a" represents the "short leg" of the 


strained plate 


ae fe a>b 
C-a; 
amb 
[cs 
end 


* Calculate deflection at which plate falls, per Van Mates 


extension to Jones 
der” 16m="20. 452%) 
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MPeser derlection counter to zero tor. this mew membrane: 
a O 


o 


$%%%%S%% Begin time-step routine %%%%%%%% 
while abs(defl) « abs(defl lim) 
if relvel « endvel 
break 
end 
% Calculate new postions and rotations at end of time step 


time = time + step; 

iS LAA «step: 
LOS YL + TIO) Step, 
X2 x2 c ZA O 
HS O AS tee 


omegal = omegal + omega _dotl*step; 
omega2 = omega2 + omega dot2*step; 
M T SC 
Pen = Pen + dPen; 
if Pen > Pmax 
Pmax = Pen; 
end 
Lek, > LBP 
f= oes 
break 
end 
IO 
L=0; 
break 
end 


if Pen » Beam 
Pen - Beam; 
break 
end 
if Pen < 0 
break 
end 


* Calculate relative translation, penetration, and change in 
impact point in this time step 

51 = [TI(1)-((BBP/2)-L)*omega doti*sin(omecgal). 
T1(2)+((LBP/2) -L) *omega_doti*cos (omegal) J ; 


S2 - (T2(1}) - (LBP2/2})*omega dot2*sin(omega2), 
T2 (2) +(LBP2/2) *omega_dot2*cos (omega2) ] ; 
reltrans"="(S2-Sil)) 4sten, 


3 Calculate direction or relativa iransla Ion 
Zetas="aban2 (re lemans (2) relervans (an 


* Calculate penetration and change in location 


aPen = sqrt (( rel trans (tea 
(reltrans(2))^2)*cos(omegal + (3*pi/2) - zeta); 

dL - sqrt((reltrans(1))^2 + (reltrans(2))°2)*sin(omegal 
EA - zeta); 
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$ Calculate the membrane deflection 
defl = defl + reltrans (2) +((LBP/2) -L) *sin(omegal) ; 


2 Calculate the resistance force GE elemen rene 


if defl < 0 
Emem - sigma y*t plate*breadth*defl^2/spacing; 
Fmem = Emem/abs (defl); 
else 
break 
end 
1£ defl < de=l lim 
deri der RE 


Emem - sigma y*t plate*breadth*defl^2/spacing; 

Fmem = Emem/abs (defl); 

EUDEUDe SA pN 

end 

% Calculate the force resulting from the "Minorsky 
mechanism" 
Rtt is the total "resistance factor" 
dRt is the differential resistance factor for this 
time step 
depth is the distance of penetration 
ddepth is the differential distance of penetration 
t 1s the aggregate in-plane structure thickness 
alpha is the bow half-entrance angle 
This formula accounts for the triangular bow wedge 


o? o? 


geometry 
and dynamic collision angle, but places no 
limits on striking ship beam. Should be modified so 
that if width exceeds beam, 


ov? H- oy? oe (T oe? o? o? ov oe 


remaining area is rectangular.... 


ddepth = sqrt (reltrans(1) 2 = reltrans (2275 

Rtt = (depth 2) *t*tan(alpha*pi/180) /(1- 
((tan(alpha*pi/180))^2/((tan((omegal-omega2)*pi/180)^2)))); 

depth = depth + ddepth; 

dRt = (depth*2)*t*tan(alpha*pi/180) /(1- 
((tan(alpha*pi/180))^2/((tan((omegal-omega2)*pi/180)^2))))- 
abs(dL)*Pen - Rtt; 


* Calculate the delta-KE from the Minorsky relation; 
Emin - dKE*10 6*dRt; 


* The corresponding force is 


o 


Fmin - Emin/abs (ddepth); 


Calculate resulting accelerations from the membrane force 
For Phase I, the membrane force is assumed to act 
perpendicularly to the hull surface 

of the struck ship, and the Minorsky force acts to oppose 
the direction of relative motion. 


o? o? 


o 


Sie 1 (Struc ent) 
Calculate the virtual mass 
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o? oe 


% Anet calculate the angle of resultant T Orce 
compared to the ship principal axis 


zetal = pi/2; % the membrane force is 
oNMpays normal to the Struck ship..... 
% Now, calculate the added mass matrix based on this 
angle 

Al = .mi*((all*einizetall 2. a22t60s (zeta 2), 
((all-a22)*sin(zetal) *cos(zetal)); ((aii- 
a22 1 "Sin (zetalp*ecos (zeta (ari*sintzetal) 2 + 
a2o*cosizetal 2); 
% Combine with physical mass to get the virtual mass 
matrix 

VM1 = M1 + Al; 
% The acceleration of translation in the global X 


Coordinate i3.: 
aXlmem = Fmem*sin(omegal) /VM1 (1,1); 


% The acceleration of translacion 1n tiend OBa: 
COOldinate as: 
aYlmem = -Fmem*cos (omegal) /VM1(2,2); 
% Calculate the angular acceleration about the ship c.g.: 
$ the cuorrenticenetace point TEE 
CP [X2 - ((LBP/2)*cos(omega2)), ((LBP/2)- 


L)*sin(omega1)]; 
% so the arm that the force acts through is: 
armi=esgre (X1 -CPP 2 RO E O E 
4 SO tChevangular Acceleration, is: 
omega dotdotimem = -((0.5-Lnd) /abs(0.5- 
Lnd))*Fmem*arm/J66V1; 


* Now calculate accelerations due to the Minorsky 
interaction, 
+ The Minorsky force is assumed to act in the direction 


opposite of relative motion. 
+ Since this force is Ine il erene ae ee Muss 


recalculate the virtual mass 


% First, calculate the angle of resultant force 
compared to the Ship principal axils 

zetal = zeta - omegal; + the angle of 
Minorsky force to thewstruchabm p M 
% Now, calculate the added mass matrix based on this 
angle 

Al = ml*{(all*sin(zetal) “2 + a22*cos(zetal) 2), 
((all-a22)*szn(zetati)*coslzeta (ace. 
aZz2)*sSin(zetal)*cos (Zebal ae a oe eee 
azz*cos(Zetad) 2) |; 
% Combine with physical mass to get the virtual mass 
matrix 


VM1 = M1 + Al; 
% Now calculate the accelerations due to the Minorsky force: 
aximin =- BEmurne*eces cd AME 
aYlmin = Fmin*sin(zeta) /VM1 (2,2) ; 
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omegal dotdotmin - -((LBP/2)-L)*Fmin*sin(zeta - pi - 
omegal)/J66V1; 


Sum the accelerations due to membrane and Minorsky for the 


toral acclera Ton due 
* due to relative motion and interaction for this time step 


axi = axlmem + axlimin; 

aYl = aYlmem + aYlmin; 

omegal dotdot = omega dotactimenm + omega lMactdo cantas 
SI DE 


o9 ov 


Calculate the virtual mass. for ene mE 
acceleration: 


% First, calculate the angle of resultant force 
compared to the ship principal axis 

zeta2 = omega2 - (omegal + pi/2) 
% Now, calculate the added mass matrix based on this 
angle 

A2 = m2*{(all*sin(zeta2) "2 + a22*cos[zeta2) 209 
((all-a22) *sin(zeta2) *cos(zeta2)); ((alil- 
a22)*sin(zeta2)*cos(zeta2)), (all*sin(zeta2)"2 + 


a22*cos 0590 m2) [I 


% Combine with physical mass to get the virtual mass 


matris 
Vile = M2 + E 


% The acceleration of translation in the global X 
coordinate is: 
aX2mem = -Fmem*sin(omegal) /VM2(1,1); 
% The acceleration of translation in the global Y 
coordinate is: 
aY2mem = Fmem*cos (omegal) /VM2 (2,2) ; 
5 The angular acceleration about the ship c.g. 1s: 
D OT dotdotmem = (-Fmem*sin(omegal-omega2 - 


pi/2) * (LBP2/2) )/J66V2; 


% Now calculate accelerations due to the Minorsky 


interaction; 
o , a * » , 
* The Minorsky force is assumed to act in the direction 


opposite of relative motion. 
* Since this force is in a different direction we must 


recalculate the virtual mass 


6 First, calculate the angle "Gr resultant iene] 
compared to the ship principal axis 

zeta2 = pi - omega2 + zeta; * the angle of 
Minorsky torce/ to he sts suo M 
Z Now, calculate the added mass matrix based on this 
angle 

A2 = m2*[(all1*sinízeta2) "2 + a02*cos (Zeta. 
((all-a22)*sin(zeta2)*cos(zeta2)); ((all- 
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ag?) *simizetaZ)*cos(zeta2)), (aàii*sin(zeta2)^2 + 
a22*cos(zeta2)~2)]; 
% Combine with physical mass to get the virtual mass 
matum x 
VM2 = M2 + A2; 

% Now calculate the accelerations due to the Minorsky force: 

aX2min = Fmin*cos (zeta-pi)/VM2(1,1); 

aY2min = Fmin*sin(zeta-pi) /VM2 (2,2); 

omega2 dotdotmin = (LBP2/2)*Fmin*sin (omega2- 
zeta) /J66V2; 


% Sum the accelerations due to membrane and Minorsky for the 
total accleration due 
* due to relative motion and interaction for this time step 


aX2 = aX2mem + aX2min; 

aY2 - aY2mem + aY2min; 

omega2 dotdot = omega2_dotdotmem + omega2 dotdotmin; 
% Calculate new velocities 

Tl uz Tuc a stepe 


qu EET 


(2) + aYi*step; 
To qo 
( 


1) + aX2*step; 
T2(2) s T2(2) + ay2*seco, 
omega doti - omega doti + omegal dotdot*step; 
omega dot2 = omega dot2 + omega2 dotdot*step; 


relvel = sqre(reltrenet i) 20 re liter I O 
Babs = Eabs + Emin + Emem - Emem last; 
Emem last = Emem; 


% Determine which tanks have been breached 
doublecargo 


end 


% double6.m 


* Script to perform time-domain analysis for double-hullea 


tanker collision: 
+ The "6" indicates that this script is Lorm hose n E e 


collision. whieh ws 
%+ from the time the inner cargo bulkhead ruptures until the 


collisrton tenas: 


Input: all the dynamic variables from double5.m 


% dt from write.m 
$ MVI from write.m 
% V2 from write m 
% VM1 from 
% VM2 thom 
% alpha from write.m 


% Output: Generated variates and values 


% Date created: 3/10/98 
* Last updated: 4/23/98 


6%%6%%S%% Begin time-step routine %%3%%%%%%% 
while relvel « endvel 
Fabs = Eabs + Emin; 


% Calculate new postions and rotations at end of time step 
time = time + step; 
A AS eos 
Yi = (MIST) stem. 
A E Dr Step, 
Y2 = 1206 TADA Steps 
omegal = omegal + omega dotl*step; 
OmegaZ = (em -Gga2 omega dotz*sutep 
[iy =a, ie 
Pen = Pen + dPen; 
if Pen > Pmax 
Pmax = Pen; 
end 
lity ee © 
Ie 
break 
end 
LE LISAS 
LS GBE; 
break 
end 
if Pen > Beam 
Pen = Beam; 


break 
end 
if Pen « O 
break 
end 


$ Calculate relative translation, penetration, and change in 
impact poin: in Chis time step 

= S1 is thetatal velocity (from linear and rotational 
mocioni of che Smpoact Pein meN Ship I. 

* S2 is the same velocity for Ship 2. 


S1 = [(T1(1)+((LBP/2) -L) *omega_dot1*sin(omegal) , 
T1(2)+((LBP/2) -L) *omega_dotl1*cos (omegal) J ; 


S2 - [T2(1) - (LBP2/2)*omega dot2*sin(omega2), 
T2(2)- (LBP2/2) *omega dot2*cos(omega2)]; 
reltrans rg estem 


* Calculate direction of relative translation 
Zeta = atan2(reltrans (2), ,reltrans(1)); 


* Calculate penetration and change in location 


dPen - sqrt((reltrans(1))^2 + 
(reltrans(2))^2)*cos(omegal « (3*pi/2) - zeta); 

dL - sqrt((reltrans(1))^2 « (reltrans(2)) 2)*sin(omega.l 
3001/2) oe tal. 


% Calculate the force resulting from the "Minorsky 
mechanism" 

Rtt as the total Iresi shance iae or" 

dRt is the differential resistance factor for this 
time step 

depth is the distance of penetration 

ddepth is the differential distance of penetration 
t is the aggregate in-plane structure thickness 
alpha is the bow half-entrance angle 

This formula accounts for the triangular bow wedge 


o? o? 


o? o9 o9 o9 oe 


Seeme try 
% and dynamic collision angle, buc pius suo 
% limits on strikiñg ship beau. Titoes ucin eena 


be modified so that if width exceeds beam, 
* remaining area is rectangular.... 


ddepth = sqrt (reltrans(1) 2 * relt*tans(2) 29€ 

Rtt - (depth^2)*t*tan(alpha*pi/180)/(1i- 
((tan(alpha*pi/180))^2/((tan((omegal-omega2)*pi/180)^2)))); 

depth = depth + ddepth; 

dRt =  (depth^2)*t*tan(alpha*pi/180)/(1- 
((tan(alpha*pi/180))^2/((tan((omegai-omega2)*pi/180)^2))))-* 
abs (dL) *Pen - Rtt; 


* Calculate the delta-KE from the Minorsky relation; 
Ena ESO Sa 


% The corresponding force is 
Fmin = Emin/abs(ddepth) ; 


Sooo Sta sHIp) 


% The acceleration of translation from membrane force 11 
the global X coordinate is: 

aX1mem = 0; 
% The acceleration of translation from membrane force in 
the global Y coordinate is: 

avYinens= 720e 


Calculate the angular acceleration about the ship c.g.: 
+ the current cota SE 
CP = [X2 - ((LBP/2)*cos(omega2)), ((LBP/2)- 
L) *sin(omegal) J ; 
* so the arm that the force acts through is: 
arm = sqrt ((XI-CPID)] 2 x (Y ICE 2a 
* so the angular acceleration due to membrane force is: 
omega dotdotıimem = 0; 


o? 


* Calculate accelerations due to the Minorsky interaction; 
% The Minorsky force is assumed to act in the direction 
opposite of relative motion. 

% Since this force is in a different direction we must 
recalculate the virtual mass 


% First, calculate the angle of resultant force 
compared to the ship principal axis 

zetal = zeta - omegal; * the angle of 
Minorsky force to tie struck Sup 
% Now, calculate the added mass matrix based on this 
angle 

Al = mil*{(all*sin(zétal) 29+ a22*cos (zen ai 
((al1-a22) *sinízetal)*cos(zetal)); ((al1- 
a22)*sin(zetal) *cos(zetal)), (all*sin(zgetal) 2 + 
a22*cos(zetal)^2)]; 
% Combine with physical mass to get the virtual mass 
matrix 


MN UMS Al; 
* Now calculate the accelerations due to the Minorsky force: 
aXimin = Fmin*cos (zeta) /VM1 (1,1); 
aYlmin = Fmin*sin(zeta) /VM1(2,2); 
omegal dotdotmin = -((LBP/2) -L) *Fmin*sinizeta i E 
omegal)/J66V1; 


* Sum the accelerations due to membrane and Minorsky for the 


total accleration due 
% due to relative motion and interaction for this time step 


axl = Yaximem + aximin-: 

aYl = aYimem + aYimin; 

omegal dotdot = omega _dotdotlmem + omegal_dotdotmin; 
< Ship 2 
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% RMismeceelerartiton of translation due to membrane Eorce 


O 


in the global X coordinate is: 


ax2mem = 0, 
a The acceleration of translation due to membrane force 
in the global Y coordinate is: 

aY2mem = 0; 
% The angular acceleration about the ship c.g. due to 
membrane force is: 

omega2 dotcotmem = 0; 


Now calculate accelerations due to the Minorsky 


interactions 
The Minorsky force is assumed to act in the direction 


opposite of relative motion. 
Since this Torce Tsai a dillere di rection MU Sis 


pecalculate the Virtue” mass 


o? H- o 


% First, calculate the angle of resultant force 
compared to the Ship principal axis 

zeta2 = pi - omega2 + zeta; % the angle of 
Minorsky Orce 10 Cho See iia some ae 
% Now, calculate the added mass matrix based on this 
angle 

A [a li*sint zeta?) 24a 2 * eos Zea 0 ey 
((a11-a22)*sin(zeta2)*cos(í(zeta2)); ((aii- 
a22)*siní(zeta2)*cosí(zeta2)), (ali*sin(í(zeta2) 2 + 
a22*cos(zeta2) ^2)]; 
% Combine with physical mass to get the virtual mass 
matrix 


VM2 = M2 + AZ; 
* Now calculate the accelerations due to the Minorsky force: 
aX2min = Fmin*cos(zeta-pi) /VM2 (1,1) ; 
ay2min = Fmn* sand zepe-pl a NS 
omega2 dotdotmin = (LBP2/2) *rmin*sin (omega2- 
zeta) /J66V2; 


* Sum the accelerations due to membrane and Minorsky for the 


total accleration due 
* due to relative motion and interaction for this time step 


ax2 = .axX2mem + ex2min- 

aY2 = aY2mem + aY2min; 

omega2 dotdot - omega2 dotdotmem « omega2 dotdotmin; 
$ Calculate new velocities 

TIAS ees ) + aXl*step; 

LIS TL) e a Yi ske. 

E 12 (2 sear 

T2 2) ccu DM REID 


necs dea = omega dO ome oo One 
omega_dot2 = omega_dot2 + omega2 _dotdot*step; 
relvel = sgqrt(reltrans(1)72 + reltrans(2)72)/step; 
Eabs = Eabs + Emin + Emem - Emem last; 
Emem last - Emem; 

* Determine which tanks have been breached 
doublecargo 

end 


mf 





% — doublecargo.m 


* the script to determine which cargo bulkheads have been 


breached during 
* each time step in the simulation for the double hull 


o 


tanker 


Datreueueated-) 3 JU OS 


oV oe 


Last revision: EA SS 
* Inputs: rana Len ron the Connon Phase «sami pire 
CUBE 


o o? 


Flags corresponding to the cargo compartments that will 
release oil 


- 1; 
while L » DBH(j) 

AP =e 
end 6 at this end,” BH(]) 1S the first bülkhecea Tar mor 
the damage location 


gc gh Sa 
oe rupture zo ae 
Tanki = sil 
if Pen >=24 


Tank1S =1; 
end 
end 
end 
ate S n 
Trop E ME 
Denc2> = Ue 
E cl ME 
Tanko = 
end 
end 
end 
IL o eea; 
O EODEM 
Tank3P = 1; 
1 -euipture 3 Gs 50 
Fania E EI 
end 
end 
end 
bao E 


if ruptures] 
Tank CAS 


end 
au 


end 
ENS TO 


end 
LRI 


end 


enc 


SS ($5 
it FUPprturez > 
Tank5P = 1; 
IE ua 
Tankse 
end 
end 


==: e 
lf rspture2 0 
DISP =I; 
If rupture) 
Tantas e 
end 
end 


== MOS 
Alpe > M0 
Slopik =r; 
MER 
Tanx6C 
end 
end 
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% iotd1.m 


% Script to perform time-domain analysis for intermediate 


o 


Oll=tignt decli tanker collision. 
$ The "i" indicates that this script is for phase 1 of the 


eollisrorn mwhrten 18 
% from the time of impact until the shell membrane ruptures. 


o 


c lnput.: aie from energy.m 
% ul from energy.m 
% V2 from energy.m 
% VM1 from energy.m 
Z VM2 from energy.m 
% alpha from vargen.m 


Output: Generated variates and values 


o 


Date created: 3/2/98 
Last updated: 4/23/98 


o? o? 


oe 


Reset flags for tank breaching to zero: 
TankiP - 
Tankis = 
TankiBOT 
Tank2P = 
tank 2 5. = 
Tank2BOT 
Man po 
Tank3S = 
Tank3BOT 
Tank4P = 
Tames ms 
Tank4BOT 
Pan 
Tank S = 
Tank5BOT 
Tanker = 
Tanke S = 
Tank6BOT 
Slop hi = 


"os ^" e ~. we o 
O e 
me e 


EM coca qp «cse 6965 NEO 
© O © © 


e 


* Determine nearest transverse structures, and distance to 


each for use in applying 
* Van Mater's extension to Jones method... 


J= i; 
while BAJ) = L 

EE 
end 
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end 
% Ensure variable "a" represents the "short leg" of the 


strained plate 


DE a»b 
C=a; 
EIE 
bzc; 
end 


% Calculate deflection at which plate faris, per Van Mate 


extension to Jones 
defl _lim= 092520 


+ Initialize new variables (Subsemipe il —) struc Sao” 
subscribe 

time = O0; 

Deme 

iE 


X2- (LBP/2)-L-*(LBP2/2)*cos (phi*pi/180); 
Y2= (Beam/2)+(LBP2/2)*sin(phi*pi/180); 
omegal=0; 

omega doc So, 
omega2=(phi+180)*(pi/180); 

eme ga JOC 20 


defl = 0; 

TL = vell, 

JUPE vel 
Fupeure e 05 
EMO Mc 
Indre e 
EIE 

eimi meo 
SAO 
A O 
Emin = 0; 
Eabs = Q; 
Emen ase anor 
Gin =) Or 

Pen = Q; 

dPen = 0; 
Pmax =a 


o? 


233605 Begin time-step routine %%%%%%%% 
while abs(detfl) < abs(defl lim) 
* Calculate new postions and rotations at end of time step 
time = time + step; 
Xa s X3 +271) step. 
MU YI ow TIH?) Sron, 
Xo c- X2 qD20 "step. 
MoNÉUY2. sp 1212) *seep-, 
omegal = omegal -« omega dotl*step; 
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omega2 = omega2 + omega_dot2*step; 
Item SIE: 

Pen = Pen + dPen; 
if Pen » Pmax 
Pmax = Pen; 

end 
if L > LBP 
L=LBP; 
break 
end 
MEE T, mO 
lE 
break 
end 
if Pen » Beam 
Pen - Beam; 
break 
end 
if Pen « O 
break 
end 


*$ Calculate relative translation, penetration, and change in 
impact point in this time step 

Si - [LEO X(ibBbb/2)Iusomedgawlotl^sumpmesc 
Ti (2)+((LBP/2) -L) *omega_dot1*cos (omegal) ] ; 


S2 = [T2(1) - (LBP2/2)*omega dot2*sin(omega2), 
T2 (2) + (LBP2/2) *omega_dot2*cos (omega2) ] ; 
rel:rans = (S2 MSIE 


* Calculate- -direction ot relative e ransler ren 
Zeta = atan2 (relitranet2) relrrans ne 


$ Calculate penetration and change in location 
dPen - sqrt((reltrans(1))^2 + 


(reltrans(2))^2)*cos(omegal + (3*pi/2) - zeta); 
dL = sqrt ((reltrans({(1)) 2 + treleranci2)) UAI NS 
+ (3*pi/2) - zeta); 


% Calculate the membrane deflection 
defl = defl + reltrans(2)+((LBP/2) -L) *sin(omegal) ; 


$ Calculate the resistance force of the membrane 
if defl « 
Emem - sigma y*t plate*breadth*defl^2/spacing; 

Fmem Emem/abs (defl) ; 


Ə 


else 
break 
end 
if defl « defl lim 
dep cem in, 
Emem - sigma y*t plate*breadth*defl^2/spacing; 
Fmem = Emem/abs (defl); 
rupture Ser 
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end 


% Calculate the force resulting frémy ene! nous, 
mechanism" 

Rett is the total “resistance ia Coi 

dRE is the differential resistance factom for ears 
tep 

depth is the distance of penetration 

ddepth is the differential distance of penetration 

t is the aggregate in-plane structure thickness 

alpha is the bow half-entrance angle 

This formula accounts for the triangular bow wedge 
geometry 
and dynamic collision angle, bue placesune 
limits on striking ship beam. Should be modified e 
that if width exceeds beam, 

% remaining area is rectangular.... 


o? o? 


time 


o? o? oV oV oV [f 


o? o? 


ddepth - sqrc (reltrans(1)^2 -« reltrans(2)^2); 

Rtt - (depth^2)*t*tan(alpha*pi/180)/ (1- 
((tan(alpha*pi/180))^2/((tan((omegai-omega2)*pi/180)^2)))); 

depth = depth + ddepth; 

CI ep cR EN. E 
((tan(alpha*pi/180))^2/((tan((omegal-omega2)*pi/180)^2))))-« 
abs (dL) *Pen - Rtt; 


% Calculate the delta-KE from the Minorsky relation; 
Emin = QKE*10°6*dRt; 


* The corresponding force is 
Fmin = Emin/abs (ddepth) ; 


$ Calculate resulting accelerations from the membrane force 
% For Phase I, the membrane force is assumed to act 
erpendicularly to the hull surface 

of the struck ship, and the Minorsky force acts to oppose 
the direction olrelative moCion. 


o 


EM Ste SE) 
Calculate the virtual mass 
First, calculate the angle of resultant force 


compared to the ship principal axis 


o? o? o? 


AA % the membrane force is 
always normal torehe stunwiek shite me 
% Now, calculate the added mass matrix based on this 
angle 

Al s mi*[(all*sin(zetat) 20r a29*cos ze A 
((all-a22)*sin(zetal) *cos(zetal)); ((ali- 
a22)*sin(zetal)*cos(zetal)), (all*sin(zetal) 2 + 
a22* COs (zeta) a2) | 
% Combine with physical mass to get the virtual mass 
Mat Y 1 


VM1 = M1 + Al; 
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% The acceleration of translation in the global X 


coordinate is: 
aXlmem = Fmem*Sin(omegal) /VM1 (1,1) ; 


% The acceleration of translation in the global Y 
coordinate is: 
aYlmem = -Fmem*cos (omegal) /VM1 (2,2) ; 
% GalculatcesunegxduasucmEccceleratmeon 2595 EOS E c 
s ehe current COn act pOlnt iS: 
CP = [X2 - ((LBP/2)*cos(omega2)), ((LBP/2) - 


L) *sin(omegal) ]; 
* SO the arm that the force acts through is: 
cum rE | (Xi eP (i A CDS 
% so the angular acceleration is: 
omega dotdotimem - -((0.5-Lnd)/abs(0.5- 
Lnd))*Fmem*arm/J66V1; 


* Now calculate accelerations due to the Minorsky 


interaction: 
* The Minorsky force is assumed to act in the direction 


Oo 
Opposite of relative motion. 
% Since this force is ina differen dlurecrtion we mus 


recalculate the virtual mass 


% First, calculate the angle of resultant force 
compared tomes Sharon) ninel loca to 

peta let o Me < the angle e: 
Minorsk tores to the secuela 
5 Now, calculate the added mass matrix based on this 
angle 

Al - ml*[(all*sin(zeta1)^2 + a22*cos (zetal)”"2), 
((all-a22) *sin(zetal) *cos(zetal)); ((all- 
a22)*szn(zetal)*cos(zetcal oUm] E? 
a22*cos(zeta1)^2)]; 
% Combine with physical mass to get the virtual mass 
matrix 


VM1 = M1 + Al; 
* Now calculate the accelerations due to the Minorsky force: 
axXimin =" hMin*cest Zeta. VMI POr 
aYlmin = Fmin*sin(zeta) /VM1(2,2) ; 
omegal dotdotmin - -((LBP/2)-L)*Fmin*sin(zeta - pi - 
omegal) /J66V1; 


* Sum the accelerations due to membrane and Minorsky for the 


total accleration due 
* due to relative motron and interact to tov sents gemma ws tcp 


axl = aXlmem + aXlmin; 

aYl = aYlmem + aYlmin; 

omegal dotdot = omega dotdotlmem + omegal dotdotmin; 
9. 9. 9- 9-9. 9. 9.9. 9.9.9.9. 9. 9.9.9. 9.9. 9. 9.9.9.9. 9. 9. 99-9. 9.9. 9. 9.0.9.9. 9.9. 99.9. 9.9.9.9. 9.9.9. 9. 9. 9. 9. 9. 
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Calculate the virtual mass for the membrane force 
acceleration: 
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% First, calculate the angle of resultant force 


compared to the ship principal axis 
zeta2 = omega2 - (omegal + pi/2); 


% Now, calculate the added mass matrix based on this 
angle 

A2 = m2*[(all*sin(zeta2)"2 + a22*cos(zeta2) 2), 
(({all-a22) *sin(zeta2) *cos(zeta2)); ((all- 
a22)*sin(zeta2) *cos (zeta2)), (all *sim(geeaz 2 
aZz2*cos(geral) 12.0% 


$ Combine with physical mass to get the virtual mass 


matrix 
VM2 = M2 + A2; 


% The acceleration of translation in the global X 
coordinate is: 

aX2mem = -Fmem*sin(omegal)/VM2 (1,1); 
$ The acceleration of translation in the global Y 


coordinates 
aY2mem - Fmem*cos (omega1)/VM2(2,2); 


% The angular acceleration about the ship c.g. is: 
omega2 dotdotmem = (-Fmem*sin(omegal-omega2- 
pi/2)*(LBP2/2))/J366V2; 


% Now calculate accelerations due to the Minorsky 


interaction; 
Qo * , * * . 
* The Minorsky force is assumed to act in the direction 


Opposite of relative motion. 
$ Since this force 1s in a different direction we must 


recalculate the virtual mass 


$ First, calculate the angle of resultant force 
compared to the ship principal axis 

zeta2 - pi - omega2 + zeta; + the angle of 
Minorsky force fo “her ah Cing onna 
% Now, calculate the added mass matrix based on this 
angle 

A2 2 m2*[l(all*sin(zeta2) 2 4 a22*cos'lzet- m S 
((all-a22) *sin(zeta2)*cos(zeta2)); ((all- 
a22)*siní(zeta2)*cos(zeta2)), (alii*sin(zeta2) 2 + 
a22*c0s (zeta) 271; 
$ Combine with physical mass to get the virtual mass 
matrix 


MIT EIN MIA 
% Now calculate the accelerations due to the Minorsky force: 
aX2min = Fmin*cos(zeta-pi) /VM2 (1,1); 
aY2min = Fmin*sin(zeta-pi) /VM2(2,2); 
omega2 dotdotmin = (LBP2/2) *Fmin*sin (omega2- 
zeta) /J66V2; 


% Sum the accelerations due to membrane and Minorsky for the 


total accleration due 
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% due to relative motion and interaction for this time step 


ov 


end 


aX2 = aX2mem + aX2min; 
aY2 = aY2mem + aY2min; 
omega2 dotdot = omega2_ dotdotmem + omega2_ dotdotmin; 


Calculate new velocities 
Sea = E + aXl*step; 


T2 Y SEO 
T2 ur smod ax2*step: 
) + aY2*step; 


Wo2 = 5 
omega doti - omega doti + omegal_ dotdot*step; 
omega dot2 - omega dot2 + omega2 p 


Loo] c E NE cuc: 
Habs = Fabs + Emin + Emem - Emem last; 
Emem lascw- Emem; 


Determine which tanks have been breached 
iotdcargo 
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% 10td2.m 


% Script to perform time-domain analysis for IOTD tanker 


Gels some 
* The "2" indicates that this script is for phase 2 of tche 


Gell perenne whieh IS 
% from the time the outer shell membrane ruptures until the 


inner bulkhead 4 s¥ecemeacecar 


Input: all the dynamic variables from iotdl.m 


% Gls from write.m 
% V1 from write.m 
% V2 From Weeen ml 
% VM1 from 
% VM2 from 
% alpha from write em 


o? 


Output: Generated variates and values 


Date created: 3/10/98 
Last updated: 4/23/98 


o? o? 


$%6%%%%% Begin time-step routine %%%%%%%%% 
while Pen < IOTDS1 

it relvel < endvel 
break 


* Calculate new postions and rotations at end of time step 
time = time + step; 

= Oe TG) Seep, 

Y= Yee Ti 2) -step- 

Ke? = D A SEE, 

¥2 = Y2 4°92 (2) seep: 


omegal = omegal + omega_dotl*step; 
omega2 = omega2 + omega _dot2*step; 
¡E 
Pen = Pen + dPen; 
if Pen » Pmax 
Pmax - Pen; 
end 
O 
Iss 
break 
end 
if L > LBP 
L=LBP; 
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break 

end 

if Pen > Beam 
Pen = Beam; 


break 
end 
if Pen < 0 
break 
end 


% Calculate relative translation, penetration, and change in 
impact point in this time step 

%4 S1 is the total veloeity Wsrom-lineer and rotational 
motion or themi mpact Point on Nio i 

% S2 is the same velocity for Ship 2. 


Sl = [T1(1)+((LBP/2) -L) *omega_dot1l*sin(omegal), 
T1 (2) +((LBP/2) -L) rre n ee 

E - [T2(1) - (LBP2/2)*omega dot2*sin(omega2), 
T2 (2) + (LBP2/2) *omega_dot2*cos (omega2) ] ; 

reltrans = (82-81) *step. 


Zeca lculate dire Cion Ol relativo trans lamen 
zeta = atan2(reltrans(2),reltrans(1)); 


% Calculate penetration and change in location 

dPen = sqrt((reltrans(1))"2 + 
(reltrans(2))0 ^2)*cosdtomedad 3b au ret 

db - sqgrt((reltranstr)) 2—* (reltrans(2))" 2) *sini(enegar 
+ (3*pi/2) - zeta); 


5 Calculate the force resulting from the "Minorsky 
mechanism" 

Rtt is the total "resistance factor" 

dRt is the differential resistance factor for this 
time step 

depth is the distance of penetration 

ddepth is the differential distance of penetration 
t is the aggregate in-plane structure thickness 
alpha is the bow half-entrance angle 

This formula accounts for the triangular bow wedge 


o? o? 


o9 o? (f o? o9? o? o? o? 


geometry 
and dynamic collision angle, but places no 
limits on striking ship beam. In the future, could 


be modified so that if width exceeds beam, 
% remaining area is rectangular.. 


ddepth = Sgrt=(reltrans (1) 27 eel tuc O 

ROCES w a mp dE a 
((tan(alpha*pi/180))^2/((tan((omegal-omega2)*pi/180)^2)))); 

depth = depth + ddepth; 

dst wro NE 
((tan(alpha*pi/180))^2/((tan((omegal-omega2)*pi/180)^2))))-* 
abs (dL)*Pen - Rtt; 
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* Calculate the delta-KE from the Minorsky relation; 
Emin - dKE*10^6*dRt; 


% The corresponding torce 15 
Fmin = Emin/abs (ddepth) ; 


* Ship 1 (Struck ship) 


% The acceleration of translation from membrane force in 
the global X coordinate is: 

aX1mem = 0; 
% The acceleration of translation from membrane force in 
the global Y coordinate is: 

ay umen <=. Om 
% Calculate the angular acceleration about the ship c.g.: 

* the current contact point is: 
CP - [X2 - ((LBP/2)*cos(omega2)), ((LBP/2)- 


L)*sin(omegal)]; 
* So the arm that the force acts through is: 
arm = sant ((X1-CP (1) ) 20 S YI PUE 
* so the angular acceleration due to membrane force is: 
omega dotdotimem = 0; 


% Calculate accelerations due to the Minorsky interaction; 
% The Minorsky force is assumed to act in the direction 
opposite of relative motion. 

% Since this force is ina different direction we must 
recalculate the virtual mass 


% First, calculate the angle ot resultant, tere 
compared to the ship principal axis 

zetal = zeta - omegal; % the angle of 
Minorsky force tothe struck chir TER: 
% Now, calculate the added mass matrix based on this 
angle 

Al = m1l*[(a@ll*sin(zetal)”"2 4+ a2Z2*cos(zetal) 2 
((all-a22)*sin(zetal) *cos (Zetal))> (lal. 
a22)*sin(zetal)*cos(zetal)), (all*sin(zetal)^2 + 
a22*cos(zetal) 2)]; 
% Combine with physical mass to get the virtual mass 
matrix 


VM1 = M1 + Al; 
* Now calculate the accelerations due to the Minorsky force: 
aXlmin = Fmin*cos(zeta)/VM1(1,1); 
aYlmin = Fmin*sin(zeta) /VM1 (2,2); 
omegal dotdotmin = =((UBP/2)-L) *Pmin*sin(a eta )] ie 
omegal) /J66V1; 


% Sum the accelerations due to membrane and Minorsky for the 


total accleration due 
% due to relative motion and interaction for this time step 


axl = aXlmem + aXlmin; 
170 


aYl = aYlmem + aYlmin; 
omegal dotdot = omega _dotdotlmem + omegal _dotdotmin; 


6 SHE D 2 
% The acceleration of translation due to membrane force 
ijeene "lobal SX eeorcinace de 

ax2mem = 0; 
% The acceleration of translation due to membrane force 
in the global Y coordinate is: 

aY2mem = 0; 
% The angular acceleration about the ship c.g. due to 
membrane force is: 

omega2 dotadorm-m mor 


% Now calculate accelerations due to the Minorsky 


Bnteractiorny 
* The Minorsky force 1S assumed to act in the direction 


opposite of relatives motion. 
2 Since this for e iS in a dirrerent Oire cCtCion Cms» 


recalculate the virtual mass 


% First, calculate the angle of resultant force 
compared to the ship principal axis 

zeta2 = pl - omega2 + zeta; * the angle of 
Minorsky tores GeO. tile Sex iain chum: m 
= Now, calculate the added mass matrix based on this 
angle 

A2 = m2*[(all*sin(zetťta2) 2" F a22*c0S (zcr E 
((all-a22) *sin(zeta2) *cos(zeta2)); ((all- 
a22)*sin(zeta2) *cos(zeta2)), (all*sin(zeta2) 2 + 
a22*cos(zeta2) ^2)]; 
% Combine with physical mass to get the virtual mass 
matrix 


VM2 = M2 + AZ; 
* Now calculate the accelerations due to the Minorsky force: 
aX2min = Fmin*cos (zeta-pi) /VM2 (1,1); 
aYv2min = Bits (zeta=p1)/VM2 (2725 
emega2 dotdocmin = (LBP2/2)"bmin*seaomeqas. 
zeta) /J66V2; 


+ Sum the accelerations due to membrane and Minorsky for the 


total accleration due 
* due to relative motion and interaction for this time step 


aX2 = axX2mem + aX2min; 
aY2 = aY2mem + aY2min; 
omega2 dotdot = omega2 dotdotmem + omega2_dotdotmin; 
% Calculate new velocities 
BUS AA ean Step: 
Lm - Lr D m Eo 
241) — T2 wake Seco, 
c RE n + aY2*step; 
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end 


omega dotl - omega dotl «oues ESL: 
omega dot2 - omega dot2 + omega2 dotdot*step; 


relvel = sgrtlteltrans io. felerana dM» SES 
Fabs = Eabs + Emin + Emem - Emem last; 
Emem last = Emem; 


Determine which tanks have been breached 
iotdcargo 


96 iotd3.m 


% Script to perform time-domain analysis for IOTD tanker 
collm&Enuom 

5 The "3" indweates that this Script is for paa seras trece 
collision Asirci e 

% from the time of impact on the inner shell, or cargo 
Þoundary until that 

* membrane ruptures. 


% Input: all the dynamic variables from iotd2.m 
% lic LEON ener m 
% V1 from energy.m 
% V2 Erom ene y 
$ VMi fu IDE EIS C EN 
% VM2 from energy.m 
% alpha from vargen.m 


oe 


Output: Generated variates and values 


Date created: 3/10/98 
Last updated: 4/23/98 


o? o 


yo 


s Determine nearest transverse structures on the inner 
shell, and distance to 

% each for use in applying Van Mater's extension to Jones 
met Mod... 


A 
while BH(j) < L 
DE NET 
end 
IB E 
Nr 
D= LABH] 1); 
else 
Io ETE 


end 
$ Ensure variable "a" represents the "short leg" of the 


strained plate 


Ir a>b 
C=a; 
a=b; 
ar 

end 


+ Calculate deflection at which clase yeas Mo Md E 


extension to Jones 
defl lin OA 


* Reset deflection counter to zere Ior errea memes 
defl = 0; 


o? 


6%%%%%% Begin time-step routine %%%%%%%% 
while abs(defl) < abs(defl lim) 
f relvel « endvel 
break 
end 
% Calculate new postions and rotations at end of time step 
time = time + step; 
A ALAS ce 
og DELE eee, 
O AE A A O 
VQ a0 Y ZAR ELA NS REO 


omegal = omegal + omega _dotl*step; 
omega2 = omega2 + omega dot2*step; 
CS Lees ee 
Pen = Pen + dPen; 
if Pen » Pmax 
Pmax - Pen; 
end 
IL PLS- NIEPE 
LLBP: 
break 
end 
E L A 
b= Os 
break 
end 


if Pen > Beam 
Pen = Beam; 
break 
end 
if Pen < 0 
break 
end 


6 Calculate relative translation, penetration, and change in 
impact point in this time step 

S1 - [TI(1)-((LBP/2) -L) *omega doti*sin(omegal), 
T1 (2) +((LBP/2) -L) *omega_dotl*cos (omegal)]; 


S2 - :T2(1) - (LBP2/2) *omega_dot2*sin(omega2), 
T2(2) - (LBP2/2) *omega dot2*cos (omega2)] ; 
reltrans - (S2-S1)*step; 


* Calculate direction of relative translation 
zeta - atan2(reltrans(2),reltrans(1)); 


*$ Calculate penetration and change in location 
dPen = sqrt((reltrans(1))"2 + 


(reltrans(2))^2)*cos(omegal + (3*p1/2) - zeta); 
aL = sqrt((reltrans(1))"2 + (reltrans(2)) 2)*sin (ones = 
+ (3*pi/2) - zeta); 
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% Calculate the membrane deflection 
defl = defl + reltrans(2)+((LBP/2)-L)*sin(omegal) ; 


€ Calculate the resistance force of the membrane 


LE debi xD 
Emem - sigma y*t plate*breadth*defl^2/spacing; 
Fmem = Emem/abs (defl); 
else 
break 
end 
ti defl < demi lim 
defl = defl lim; 
Emem - sigma y*t plate*breadth*def1^2/spacing; 
Fmem = Emem/abs (defl); 
BUSeUureZ = n; 
end 


% Calculate the force resulting from the "Minorsky 
mechanism" 

Retais the total "resistance tactror® 

dRt is the differential resistance factor for this 
time step 

depth is the distance of penetration 

ddepth is the differential distance of penetration 
t is the aggregate in-plane structure thickness 
alpha is the bow half-entrance angle 

This formula accounts for the triangular bow wedge 
ee cer 

and dynamic collision angle, but places no 

limits on striking ship beam. Should be modified so 
f width exceeds beam, 

remaining area is rectangular... 


o? o? o? o? o? o? o9 


o? o? 


that 


oe H- 


ddepth = sqrt (relrtransii) 2. + releeans 2) 2s 
Rtt = (depth*2) *t*tan(alpha*pi/180) /(1- 
( (tan (alpha*pi/180))^2/((tan(í((omegai-omega2)*pi/180)^2)))); 
depth = depth + ddepth; 
art = (depth 2) 4tttan(alpha-e1/ ee sere 
((tan (alpha*pi/180))*2/((tan( (omegal-omega2) *pi/180)7~2))))+ 
abs(dL)*Pen - Rtt; 


* Calculate the delta-KE from the Minorsky relation; 
Em Es oS CES 


* The corresponding force is 
Fmin = Emin/abs (ddepth) ; 


Calculate resulting accelerations from the membrane force 
For Phase I, the membrane force is assumed to act 
erpendicularly to the hull surface 
of the struck ship, and the Minorsky force acts to oppose 
the direction of relative motion. 


o? o? 


oe FO 


Ship I (Struek ship) 
Calculate the virtual mass 


p 


o o? 


% First, calculate the angle of resultant force 
compared to the ship principal axis 


zetal = pi/2; % the membrane force is 
always normal Eo the steuck Shir TSE 
% Now, calculate the added mass matrix based on this 
angle 

Al = ml*[(all*sin(zetav) 2a “cee eca as 
((all-a22)*sin(zetal)*cos(zetal)); ((aii- 
a22)*sin(zetal) *cos(zetal)), (all*sin(zetal)72 + 
a22*cos(zeta1)^2)]; 
% Combine with physical mass to get the virtual mass 
Meteo ds x 

VM1 = M1 + Al; 
$ The acceleration of translation in the global X 


Coordinates: 
aX1lmem = Fmem*sin(omegal) /VM1 (1,1); 


% The acceleration of translation in the global Y 
coordinates: 
aYlmem = -Fmem*cos (omegal) /VM1 (2,2) ; 
% Calculate the angular acceleration about the ship c.g.: 
+ PEheteurrent tesncact pestes 
CP = [X2 - ((LBP/2)*cos(omega2)), ((LBP/2) - 


L)*sin(omegal)]; 
$ so the arm that the force acts through is: 
arm = sgrt((X1-CP(1)) 2 + (YT GE I 
% so the angular acceleration is: 
omega_dotdotimem = -((0.5-Lnd) /abs(0.5- 
Lnd))*Fmem*arm/J66V1; 


% Now calculate accelerations due to the Minorsky 
interaction, 

% The Minorsky force is assumed to act in the direction 
opposite of relative motion. 

% Since this force is in a different direction we must 


o 


recalculate the virtual mass 


$ First, calculate the angle of resultantes 
compared to the ship principal axis 

fetal DG esum $ the angle of 
Minorsky force to the Sseruck siip ae 
% Now, calculate the added mass matrix based on this 
angle 

Al = mi*[(all*sin(zetal) “2 + a22*cos(zetal)~2), 
((all-a22) *sin(zetal) *cos(zétal)). ((all- 
a22)*sin(zetal)*cos(zetal)), (all*sin(zetal)~2 + 
a22*cos(zetal)~2)]; 
$ Combine with physical mass to get the virtual mass 
matrix 


VMI = Maas AL: 
* Now calculate the accelerations due to the Minorsky force: 
aXimin = Fmin*cos(zeta) /VM1(1,1); 
ayimin scBmin*sin(zeta snb 229 
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cec NEcteotuem - I UEBD/2) E) *pn Sime CI MN D 
omegal)/J66V1; 


Sum the accelerations due to membrane and Minorsky for the 


otal acoleration due 
due to relative motion and interaction for this time step 


A T o 


aXl = aXlmem + aXlmin; 

aYl = aYlmem + aYlmin; 

omegal dotdot - omega dotdotimem « omegal dotdotmin; 
SM 


o oV 


Calculate the virtual mass for the membrane force 
acceleration: 


$ First, calculate the angle of resultant force 
compared to the ship principal axis 

zeta2 = omega2 - (omegal + pi/2); 
$ Now, calculate the added mass matrix based on this 
angle 

A2 = m2*[(all*sin(zeta2)"2 + a22*cos(zeta2)"2), 
((all-a22)*sin(zeta2)*cos(zeta2)); ((all- 
a22)*sin(zeta2) *cos(zeta2)), (all*sin(zeta2)"2 - 


Epic rci NE NN 


% Combine with physical mass to get the virtual mass 


matrix 
VM2 = M2 + A2; 


% The acceleration of translation in the global X 
coordinate is: 

aX2mem = -Fmem*sin(omegal) /VM2(1,1); 
% The acceleration of translation in the global Y 


coordinate is: 
aY2mem - Fmem*cos (omegal)/VM2 (2,2); 


% They angular acceleration about. the Ship cc ME E 
clesie == menos in omega omes 
pi/2)*(LBP2/2))/J66V2; 


% Now calculate accelerations due to the Minorsky 


interaction; 
% The Minorsky force is assumed to act in the direction 


opposite of relative motion. 
= Since this forces Inoa Aa prerene ac cer lee wow iluS 


recalculate the virtual mass 


% First, calculate the angle of resultant force 
compared to the ship principal axis 

zeta2 = pi - omega2 + zeta; a tne ang leven 
MrnorskcMrorceawmbosEa ser IIS 
% Now, calculate the added mass matrix based on this 
angle 

AQ] m2* | laliesinm zetal) 2°94 29x: 7 Ze 
((all-a22) *sin(zeta2) *cos(zeta2)); ((al1l- 


ey 


a22)*sin(zeta2) *cos(zeta2)), (all*sin(zeta2)~2 + 
a22*cos (zeta2) ^2)]; 
6 Combine with physical mass to get the virtual mass 
matrix 
VM2 = M2 + A2; 
% Now calculate the accelerations due to the Minorsky force: 
aX2min = Fmin*cos(zeta-pi) /VM2 (1,1); 
aY2min = Fmin*sin(zeta-pi) /VM2 (2,2); 
omega2 dotdotmin = (LBP2/2) *Fmin*sin (omega2- 
zeta) /J66V2; 


% Sum the accelerations due to membrane and Minorsky for the 
total accleration due 
% due to relative motion and interaction for this time step 
aX2 = axX2mem + aX2min; 
aY2mem + aY2min; 


aY2 = 
omega2 dotdot = omega2 dotdotmem + omega2 dotdotmin; 


% Calculate new 
To E SIS ties) 


Tit2) v= (2) Vteay cece, 
D Cie) T2(1) + ax2*step; 
T2 (2) -= T22) ae Sree, 
omega doti - omega dotl + omegal dotdot*step; 


omega dot2 = omega dot2 + omega2 dotdot*step; 


Vvelecittes 
Prac Sew, 


relvel = sqrt (reltrans(1)"2 + reltrans (2) ~2)/sten, 
Fabs = Eabs + Emin + Emem - Emem last; 
Emem last = Emem; 

% Determine which tanks have been breached 


YO SS 


end 
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% 10td4.m 


a SCTipE tO pDErrorm cime domain analysis tos IOTDI anker 
Collision: 

* The "4" indicates that this script is for phase 4 of the 
collision, which ls 

5 from the time the inner shell membrane ruptures until the 


inner cargo bulkhead is 
* contacted. 


Input: all the dynamic variables from iotd3.m 


% Giz from write.m 

5 ND from write.m 

$ a from write.m 

% VM1 from 

% VM2 from ' 
% alpha Erom write m 


Output: Generated variates and values 


o? 


Date created: 3/10/98 
Last updated: 4/23/98 


o AS 


o? 
o? 
o 
oy? 
o? 
ov 
o? 
o? 
o 


6%%%%%S% Begin time-step routine 
hile Pen < IOTDS2 
f relvel « endvel 

break 


= 


L 


end 
Babs = Eabs + Emin; 
% Calculate new postions and rotations at end of time step 
time = time + step; 
O G ASS 
YL = Yi ree cece, 
X2 = XO ee cm 
MONDO (2 seer, 
omegal = omegal + omega _dotl*step; 
omega2 = omega2 + omega_dot2*step; 
GS G ardi 
Pen = Pen + dPen; 
if Pen > Pmax 
Pmax = Pen; 
end 
Be fee wy O 
E=0% 
break 
end 
if L > LBP 
DIES 
break 


mo 


end 
if Pen > Beam 
Pen = Beam; 


break 
end 
if Pen < 0 
break 
end 


* Calculate relative translation, penetration, and change im 
impact point in this time step 

* Si is the total velocity (from linear and rotational 
motion) of the impact point NND 

* S2 is the same velocity for Ship 2. 


Sl = (Ti (1) +((LBP/2) E omega deri sia nes ne 
T1(2)+((LBP/2) -L) *omega_doti*cos(omegal) ] ; 


S2 = (92(1) = (BBP2/ 2) omega dee27ss a emeqacy 
T2 (2) +(LBP2/2) *omega_dot2*cos (omega2) ] ; 
reltrans = (S2-S1)*step; 


$ Calculate avrection ol relative rans lacie 
zeta = atan2 (reltrans (2) ) relteransi1) e 


* Calculate penetration and change in location 
dPen = sqrt((reltrans(1i))72 + 


(reltrans(2))^2)*cos(omegal + (3*pi/2) - zeta); 
dL = sqrt ((reltrans(i))"2 + (reltrans(2))~2)*sin(omegam 
+ (3*pi/2) - zeta); 


% Calculate the force resulting from the "Minorsky 
mechanism" 

Rttwms the total "resrsrtaucewftoctoun 

dRt is the differential resistance factor for this 


o oV 


time step 

% depth is the distance of penetration 

* ddepth is the differential distance of penetration 

* t is the aggregate in-plane structure thickness 

% alpha is the bow half-entrance angle 

* This formula accounts for the triangular bow wedge 
geometry 

* and dynamic collision angle, but places no 

* limits on striking ship beam. In the future, could 


be modified so that if width exceeds beam, 
5$ remaining area is rectangular.... 


ddepth = sqrt (reltrans(2) (2°) erelrnans i UE 

Rtt = (depth^2) *t*tan (alpha*pi/180)/(1- 
((tcan(alpha*pi/180))^2/((tan((omegai-omega2)*pi/180)^2)))); 

depti = depti depu. 

GRE = (depth°2)*t*tanl(alpha Pi ICOS 
((tan(alpha*pi/180))^2/((tan((omegal-omega2)*pi/180)^2))))-« 
abs(dL)*Pen - Rtt; 
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% Calculate the delta-KE from the Minorsky relation; 
Emin - dKE*10^6*dRt; 


% The corresponding force is 
Fmin = Emin/abs (ddepth) ; 


So) Sl Op ME EEUU Smp) 


% The acceleration Of translation from membrane Torce in 
the global X coordinate is: 
aximemt= 0 


The accelera ion O -translation from momoran erce In 


ehe global Y "coord mares Lis: 
aYlmem = O0; 


Calculate the angular acceleration about the ship c.g.: 
+ rhe current eontact perme ais: 
CP = [X2 - ((LBP/2)*cos(omega2)), ((LBP/2) - 
L)*sin(omegal) J; 
% so the arm that the force acts through is: 
arm sqreQ XL-cPUI) ZUM MENU M e ME 
* so the angular acceleration due to membrane force is: 
omega docdorcimemni er, 


o? 


% Calculate accelerations due to the Minorsky interaction, 
% The Minorsky force is assumed to act in the direction 


opposite of relative motion. 
$ Since this force is in a different direction we must 


recalculate Cc PH amaS 


% First, calculate the angle of resultant force 
compared to the ship principal axis 

zetal = zeta - omegal; % the angle of 
Minorsky force to Lhe si rvek=sieoea 
$ Now, calculate the added mass matrix based on this 
angle 

A1 - mi*[(all*sin(zetal)^2 « a22*cos(zetal) 2), 
((all-a22)*sin(zetal) *cos(zetal)); ((all- 
a22)*sin(zetal)*cos(zetal)), (all*sin(zetal) 2 - 
a22*cos(zetal) 2)1; 
% Combine with physical mass to get the virtual mass 
matrix 


VM1l = M1 + Al; 
* Now calculate the accelerations due to the Minorsky force: 
aXimin = Fmin*cos (zeta) /VM1(1,1); 
aYimin = Fmin*sin(zeta) /VM1 (2,2) ; 
omegal dotdotmin = (LBP/2)-L)*Fmin*sin(zeta - pi - 
omegal) /J66V1; 


% Sum the accelerations due to membrane and Minorsky for the 
total aceleration que 
% due to relative motion and interaction for this time step 
axl = aXlmem + aXlmin; 
aYl = aYlmem + aYlmin; 
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omegal dotdot = omega_dotdotlmem + omegal dotdotmin; 


cS ng e 
% The acceleration of translation que to memerane eres 


in the global X coordinate is: 
aX2mem = O0; 


% The acceleration of translation due TComenbrane Mores 
in the global Y coordinate is: 
aY2mem*t=>+0 = 


% The angular acceleration about the ship c.g. due to 
membrane force is: 
omega2 dotdotmem = 0; 


% -Now calculate accelerations miii e vromenen)ineacc) 


interaction; 

$ The Minorsky force 1S assumed to det in che direction 
opposite of relative motion. 

$ Since this force is in a differeənt direction we muse 
recalculate the virtual mass 


% First, calculate the angle of resultant force 
compared to the ship principal axis 

zeta2 = pl - omega2 + zeta; % the angle of 
Minorsky fonce tothe striking Ship..... 
% Now, calculate the added mass matrix based on this 
angle 

A2 - m2*[(all*sin(zeta2) 2 + a22*cos({zebaZ) “2m 
((a11-a22)*sin(zeta2)*cos(zeta2)); ((a1l1- 
a22)*sin(zeta2) *cos(zeta2)), (all*sin(zeta2)7~2 + 
a22*cos (zeta2)~2)]; 
% Combine with physical mass to get the virtual mass 
matrix 


VM2 = M2 + A2; 
% Now calculate the accelerations due to the Minorsky force: 
ax2min = Fmin cosl zeta-PIP Die DE 
aY2min = Fmin*sinlzeta- pi) Ma 
omega2 dotdotmin - (DBP2/2)*Fmsn*sScne mec za 
zeta) /J66V2; 


% Sum the accelerations due to membrane and Minorsky for the 


total accleration due 
* que to relative motion and interaction for this time step 


aX2 = aX2mem « aX2min; 

aY2 = aY2mem + aY2min; 

omega2 dotdot - omega2 dotdotmem « omega2 dotdotmin; 
% Calculate new velocities 

TIAS i + aXl*step; 

T4502) RIED) exacte 

T2(1) = T2(1) + aX2*step; 

T2(2) = 7242) + <aAvYZ2*step; 


omega _ doti + omegal dotdotk gopi 
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e 


end 


omega dot2 - omega dot2 + omega2 dotdot*step; 


relvel - sqrt(reltrans(1)^2 + reltrans(2)72)/step; 
Eabs = Eabs + Emin « Emem - Emem last; 
Emem loot = Enea, 


Determine which tanks have been breached 
lotdcargo 


% 10td5.m 


* Script to perform time-domain analysis for IOTD tanker 
collision 
% The "5" indicates that this script is for phase 5 of the 


collision, whiem is 
% from the time of impact on the centerline bulkhead, until 


that 
* membrane ruptures. 


Input: Yall’ the dynamie Variables ironereiea+ aim 


$ de from energy.m 
% E LEON Fenster ml 
6 V2 from energy.m 
% VM1 ISO ES ay 
% VM2 from energy.m 
% alpha from vargen.m 


oe 


Output: Generated variates and values 


Date created: 3/10/98 
Last updated: 4/23/98 


o o 


5$ Determine nearest transverse structures on the shell, and 


distance to 
% each for use in applying Van Mater's extension to Jones 


method.... 


- 1; 
while BH(j) « L 
asi; 
end 
a = BR Li; 
LE J e 
b = L - BH(j-1); 
else 
Or D 


end 
% Ensure variable "a" represents the "short leg" of the 


strained plate 


ge a>b 
C-a; 
asp 
ccs 
end 


% Calculate deflection at which plate fails, per Van Mater's 


extension to Jones 
demi Tum 5.034 5c, 
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"Reset daenleetuon counter Co zero ftorwthis newmnembrane: 


o 


deu 07 


o 
A 
o9 
o? 
o? 
o? 
oe 
o 
o\? 


$%%%3%S%% Begin time-step routine 
while abs(defi) < abs(defl lim) 
f relvel < endvel 
break 
end 
* Calculate new postions and rotations at end of time step 
time = time + step; 
DNE SUCI UNUS SEI 
Xi =-Y1 A seep; 
XQ c X NONAS. 
pu c vcn CD 
omegal = omegal « omega dotl*step; 
omega2 = omega2 -« omega dot2*step; 
DM uS 
Pen = Pen + dPen; 
if Pen » Pmax 
Eulate 
end 
LE Lasc LBP 
L=hBe; 
break 
end 
Jof Sec. 
IS 
break 
end 
if Pen » Beam 
Pen - Beam; 
break 
end 
JN pone 0 
break 
end 
* Calculate relative translation, penetration, and change in 
impact point in this time step 
S1 = [T1(1)+((LBP/2) -L) *omega_doti*sin(omegal), 
T1(2)-((LBP/2) -L) *omega doti*cos(omegal)]; 


S2 - [T2(1) - (LBP2/2)*omega dot2*sin(omega2), 
T2 (2) + (LBP2/2) *omega_dot2*cos (omega2) ] ; 
reltrans = {(S2-Si) Stem 


6 Calculate Gixguection ot relabiv]= translation 
zeta = atan2 (reltrans(2) reltrans(1)); 


* Calculate penetration and change in location 

dPen - sqrt((reltrans(1))^2 + 
(reltrans(2))^2)*cos(omegal + (3%*pi/2) - zeta); 

dL = sqrt((reltrans (2!) ) 2 4 (reltransi2 yea emecat 
+ (3*pi/2) - zeta); 
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% Calculate the membrane deflection 
defl = defl + reltrans(2)+((LBP/2) -L) *sin(omegal) ; 


* Calculate the resistance force of the membrane 


LL Cee bread) 
Emem - sigma y*t plate*breadth*defl^2/spacing; 
Fmem = Emem/abs (defl); 
else 
break 
end 
Lf desto der mm 
defl decem. 


Emem - sigma y*t plate*breadth*defl1^2/spacing; 
Fm-m = Emem/abs (def£f1); 
ruptures NI 
end 
% Calculate the force resulting from the "Minorsky 
mechanism" 
Rtt is the total "resistance factor" 
dRt is the differential resistance factor for this 
time step 
depth is the distance of penetration 
ddepth is the differential distance of penetration 
t is the aggregate in-plane structure thickness 
alpha is the bow half-entrance angle 
This formula accounts for the triangular bow wedge 
de ometry 
and dynamic collision angle, but places no 
limits on striking ship beam. Should be modified so 
chat 1£ width exceeds beam, 
$ remaining area is rectangular.... 


oV? oV oV oV? oV? o wil va 


o? oV 


ddepth = sqrt (reltrans(1)°2 + reltrans (207m, 

Rtt - (depth^2)*t*tan(alpha*pi/180)/(1i- 
((tan(alpha*pi/180))^2/((tan((omegal-omega2)*pr/180) 2008 

depth = depth + ddepth; 

dRt =  (depth^2)*t*tan(alpha*pi/180)/ (1- 
((tan(alpha*pi/180)) “2/((tan((omegal-omega2) *pi/180) ~2))))+ 
abs(dL)*Pen - Rtt; 


* Calculate the delta-KE from the Minorsky relation; 
Emin = ae 20) 6h. 


% The corresponding force is 
Fmin = Emin/abs (ddepth) ; 


Calculate resulting accelerations from the membrane force 
For Phase I, the membrane force is assumed to act 
perpendicularly to the hull surface 

of the struck ship, and the Minorsky force acts to opr onie 
the direction of relative motion. 


o? oV 


o 


Shion X SbEU SK SES) 
Calculate the virtual mass 
186 


ov o 


% First, calculate the angle of resultant force 


o 


compared to the ship principal axis 


PORUM DE 5$ che membrane force is 
always noemi o Ene Semuck Ship..:. . 
% Now, calculate the added mass matrix based on this 
angle 

A1 - mi*[(ali*sin(zetal1) 2 4 a22*cos(zetal) 2), 
((all1-a22)*sin(zetal)*cos(zetal)); ((all- 
a22)*sin(zetal)*cos(zeta1)), (all*sin(zetal)"2 + 
a22*cos(zetal) 2)]; 
% Combine with physical mass to get the virtual mass 
matrix 

VM1 = M1 + Al; 
% The acceleration of translation in the global X 


coordinate is: 
aXimem - Fmem*sin(omegal)/VM1(1,1); 


$ The acceleration of translation in the global Y 
coordinate is: 
aYlmem = -Fmem*cos (omegal) /VM1 (2,2) ; 
% Calculate the angular acceleration about the ship c.g.: 
5 the current contact point is: 
CP = [X2 - ((LBP/2) *cos(omega2)), ((LBP/2)- 


L) *sin(omegal) J ; 
% so the arm that the force acts through is: 
Arm, = equip oc DE 
* So the angular acceleration is: 
omega dotdotimem - -((0.5-Lnd)/abs(0.5- 
Lnd))*Fmem*arm/J66V1; 


* Now calculate accelerations due to the Minorsky 


interaction; 

% The Minorsky force is assumed to act in the direction 
opposite of relative motion. 

5$ Since this force is in a different direction we must 
recalculate the virtual mass 


$ First, calculate the angle of resultant force 
compared to the ship principal axis 

zetal = zeta - omegal; 5 the angle of 
Minorsky force Bothe Struck ShI Ta un 
% Now, calculate the added mass matrix based on this 
angle 

A1 - ml*[(all*sin(zetal1)^2 + a22*cos(zetal) ^2), 
((al1-a22)*szm(zetal)*cos(zetal)); Comum 
a22)*sin(zetal) *cos(zetal)), (all*sin(zetal) 2 + 
a22*cos(zetal)~2)]; 
5 Combine with physical mass to get the virtual mass 
malos 


Vil =" Maia + Ad 
* Now calculate the accelerations due to the Minorsky force: 
aximin = Fmin*coaizeraly Wiincl wie 
aYlmin = Fmin*sin(zeta) /VM1 (2,2) ; 
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omegal dotdotmin = -((LBP/2)-L)*Fmin*sin(zeta - pi - 
omegal) /J66V1; 


% Sum the accelerations due to membrane and Minorsky for the 


Eetal-acelleration due 
% due to relative motion and interaction for this time step 


aXl = aXlmem + aXlmin; 

aYl = aYlmem + aYlmin; 

omegal dotdot = omega_dotdotimem + omegal dotdotmin; 
epe 


o? oe 


Calculate the virtual mass for the membrane ER 
acceleration: 


% First, calculate the angle of resultank tore] 
compared to the shipi princi a manis 

zeta2 = omega2 - (omegal + pi/2); 
% Now, calculate the added mass matrix based on this 
angle 

A2 2 m2*[(all*srn(izeta2) 2 + a22ee@os ( Zeral sae 
((all-a22)*sin(zeta2) *cos(zeta2)); ((all- 
a22)*sin{zeta2) *cos{zetazZ) 9,5 (ali*simevetaz) 2: 


a22*cos(zeta2) 2)]; 


% Combine with physical mass to get the virtual mass 


matrix 
VM2 = M2 F A2; 


% The acceleration of translation in the global X 
coordinate is: 

aX2mem = -Fmem*sin(omegal) /VM2(i,1); 
% The acceleration of translation in the global Y 


coordinate is: 
aY2mem = Fmem*cos (omegal) /VM2 (2,2) ; 


$ The angular acceleracion about the ship Graze. 
omegaz adoecoumem = {(-Fmem* sin (eomegal -enegaz- 
pi/2)*(LBP2/2))/d66V2; 


% Now calculate accelerations due to the Minorsky 


interaction; 
o * . . * . 
+ The Minorsky force is assumed to act in the direction 


opposite of relative motion. 
% -Since this force is in a dPPforentodrcocbrone S705 


recalculate the virtual mass 


$ First, calculate the angle of resultant force 
compared to the ship principal axis 

zeta2 = pl - omega2 + zeta; % the angle of 
Minorsky force to the striking ship..... 
% Now, calculate the added mass matrix based on this 
angle 


A2 s m2*[(ali*sin(zeta2) 2 £°a22*ees (Zeta. 2 
((all-a22)*sin(zeta2)*cos(zeta2)); ((aili- 
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ias Rao cos (zeta). (eall*saumigetal) 2 4 
a22 *GeszZerae, 2) 15 
6 Combine with physical mass to get the virtual mass 
matrix 
VM2 = M2 + A2; 
% Now calculate the accelerations due to the Minorsky force: 
aX2min - Fmin*cos (zeta-pi)/VM2(1,1); 
aY2min = Fmin*sir(zeta-pi) /VM2 (2,2); 
omega2_ dotdotmin = (LBP2/2)*Fmin*sin (omega2- 
zeta) /J66V2; 


% Sum the accelerations due to membrane and Minorsky for the 


total accleration due 
* due to relative motion and interaction for this time step 


aX2 = aX2mem + axX2min; 

aY2 = aY2mem + aY2min; 

omega2 dotdot = omega2_dotdotmem + omega2 dotdotmin; 
% Calculate new velocities 

E E ae sree, 

Pi 2) = Ei (2 ea AS O 

PI) SI A) A ae Secu, 

A E O O 


omega doti = Omega dotl r omnmcga M EIE 
omega: dot2 < omega dot2 +Tomcega2 ar Caor aa 


relvel = sqrt(reltrans(1)*2 -« reltrans(2)^2)/step; 
Babs = Eabs + Emin « Emem - Emem last; 
Emem last - Emem; 
% Determine which tanks have been breached 
RI NINOS 


end 
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% iotd6.m 


% Script to perform time-domain analysis for IOTD tanker 


col Lists e 


% Tne "6" indicates that this script is for phase 6 of the 


collision, wissen 5 


% from the time the centerline bulkhead ruptures until the 


collision ends. 


* Input: all the dynamic variables from doubles5. 
% dt from write.m 

% vI from wr Tecan 

% V2 from write.m 

e VMi from 

% VM2 from 

$ alpha from write.m 


% Output: Generated variates and values 


2 Date credis c ND oe 
* Last updated: 4/23/98 


6%6%%%%%S% Begin time-step routine %%%%%%%%% 
while relvel « endvel 


Fabs - Eabs + Emin; 
+ Calculate new postions and rotations at end of 
time = time + step; 


ee eee eo. 


Yl = Yl + Ti(2)*step; 
ORE Nodo ce EU cime 
Y2 = Y2os 2 (2) 45 Gep; 
omega = sees EGISSET 
omega2 = omega2 + omega dot2*step; 
a eee 
Pen = Pen + aPen; 
1f Pen > Pmax 
Pmax = Pen; 
end 
E AO 
ME 
break 
end 
Lt EN eee 
E=EBE" 
break 
end 
if Pen > Beam 
Pen = Beam; 
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time step 


break 
end 
if Pen « O 
break 
end 


% Calculate relative translation, penetration, and change in 
impact point ine e e Ime stem 
% S1 is the total velocity (from linear and rotational 


o 


WMOEtOn) (Caste Impact polar ron tea l. 
% 52 is the same velocity for Ship 2. 


S1] s [T1(1)-((LBP/2)-L)*omega doti*sin(omegal), 
Ti(2)-«((LBP/2)-L) ME m 

E - [T2(1) - (LBP2/2) *omega dot2*sin(omega2), 
T2(2) -(LBP2/2) *omega dot2*cos (omega2)]; 

reltrans = (S2-S1)*step; 


+ Calculatela:rection fot trelatiuve trans!acion 
zeta = atan2(reltrans (2), reltrans(1)); 


% Calculate penetration and change in location 

dPen = sqrt((reltrans(1))°2 + 
(reltrans(2))^2)*cos(omegal + (3*pi/2) - zeta); 

dL = sqrt((reltrans(1))°2 + (reltrans(2))72)*sin(omegal 
+. (3*p1/2) = “Zerae 


% Calculate the force resulting from the "Minorsky 
mechanism" 
* Rtt is the total "resistance factor" 
* dRt is the differential resistance factor for this 
time step 
depth is the distance of penetration 
idepth is the differential distance of penetration 
t is the aggregate in-plane structure thickness 
alpha is the bow half-entrance angle 
This formula accounts for the triangular bow wedge 
ry 
and dynamic collision angle, but places no 
limits on striking ship beam. In the future, could 
be modified so that if width exceeds beam, 
5$ remaining area is rectangular.... 


geome 


o? o9 (T o? d? de o? oe 


ddepth = sqrt (réltrans(i) 2 4 relerans (7) 22 

Rtt = dE or 
((tan(alpha*pi/180))^2/((tan((omegal-omega2)*pi/180)^2)))); 

depth = depth + ddepth; 

dRt = vow 
( (tan (alpha*pi/180))^2/((tan((omegal-omega2)*pi/180)^2))))-* 
abs(dL)*Pen - Rtt; 


% Calculate the delta-KE from the Minorsky relation; 
Emin = "dKE*10 6*dRt- 


191 


Sine corresponding emet 
Fmin = Emin/abs (ddepth) ; 


* Ship 1 (Struck ship) 


% The acceleration of translation from membrane force am 
the global X coordinates: 

aXimem = 0; 
% The acceleration of translation from membrane force in 
the global Y coordinate is: 

axdtmene- 
% Calculate the angular acceleration about the ship com 

+ the current cowWwboctapemup ec 
CP = [X2 - ((LBP/2)*cos(omega2)), ((LBP/2)- 


L) *sin(omegal)]; 
+ sO the arm that the rorce acts througmeis- 
arm =sgsgrt  (XI-CP( DD 2 (Yin DNE E 
% so the angular acceleration due to membrane force is: 
omega dotdotimem = 0; 


% Calculate accelerations due to the Minorsky interaction; 
% The Minorsky force is assumed to act in the direction 


opposite of relative motion. 
% Since this force ws 1n a different di rection neeme 


recalculate the virtual mass 


% First, calculate the angle of resultant Corte 
Compared EG Ene=sntp prin ipa ees 

zetal = zeta - omegal; * the angle of 
Minorsky korceUto xnsetstucuss sump. 0 
% Now, calculate the added mass matrix based on this 
angle 

Al = ml* | (alil*sin(zetal) 2 + a2Z2*cos (Zeta. M 
((all-a22) *sin(zetal) *cos(zetal)); ((al1l1- 
a220canzetsu coc MS cgu qe ee me 
a22%*c0s (zeta 2115 
% Combine with physical mass to get the virtual mass 
matrix 


VM1 = M1 + Al; 
* Now calculate the accelerations due to the Minorsky force: 
aXlmin = Fmin*cos(zeta)/VM1(1i,1); 
aYimin = Fmin*sin(zeta) /VM1 (2,2); 
omegal dotdotmin = -((LBP/2)-L)*Fmin*sin(zeta - pi - 
omegal) /J66V1; 


6 Sum the accelerations due to membrane and Minorsky for the 
otal accleration due 
* due to relative motion and interaction for this time step 
axi = aXimem + aXlmin; 
aYl = aYlimem + aYimin; 
omegal dotdot = omega_dotdotimem + omegal_dotdotmin; 


a 


SOS 2 
19>? 


% TMaemaceceteration Of translation due to membrane force 
in@eie Global Kk coordindie s 
ax2mem = O0; 


% The acceleration of translation due to membrane force 
in the global Y coordinate is: 

ayvy2men = 07 
% The ancular acceleration about the ship c.g. due to 
membrane force is: 

omega2 dotdotmem - 0; 


% Now calculate accelerations due to the Minorsky 


interaction; 
% The Minorsky force is assumed to act in the direction 


opposite of relative motion. 
Since CiS T Orce 15.1 e diilesene dalreceaion we alte 


mcd cule vir udal mass 


% First, calculate the angle of resultant force 
compared to the ship principal axis 

zeta2 - pi - omega2 + zeta; % the angle orf 
Minorsky Loree to tre Str ikinc slo Um 
% Now, calculate the added mass matrix based on this 
angle 

AQ = m2* | (all*sin(zeta2) 2 += ae22*cos(zeral) 2). 
((all-a22) *sin(zeta2) *cos(zeta2)); ((all- 
a22)*sin(zeta2) *cos(zeta2)), (all*sin(zeta2)~2 + 
a22 Cos (Zetal) 2a. 
% Combine with physical mass to get the virtual mass 
matrix 


VM2 = M2 + A2; 
% Now calculate the accelerations due to the Minorsky force: 
amina = Pmin*cos(geca-pi je VM (1 
Ami - Fmintsinl zeta P7 VM? eee ee 
omega2 dotdotmin = (LBP2/2) *Fmin*sin (omega2- 
zeta) /J66V2; 


+ Sum the accelerations due to membrane and Minorsky for the 


total accleration due 
% due to relative motion and interaction for this time step 


aX2 = aX2mem + axX2min; 

aY2 = aY2mem + aY2min; 

omega2 dotdot = omega2 dotdotmem + omega2_ dotdotmin; 
% Calculate new velocities 

Ti] oS a Sep: 

T4 (2) = Oa ee ae Se Sp: 

T2 (1) = T2 D ee ete: 

T22) ae Oe E AMS O 


omega doti omega lao omega dotado SE. 
omega dorm tome ado ome OPS a 
relvel - sqrt(reltrans(1)^2 « reltrans(2)^2)/step; 
Eabs - Eabs + Emin + Emem - Emem_last; 
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Emem last - Emem; 


% Determine which tanks have been breached 
Lotdcargo 


end 
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%  iotdcargo.m 


o 


breached during 
* each time step in the simulation for the mid-deck tanker 


o 


% the script to determine which cargo bulkheads have been 


Dace reae OS 


o? oe 


Last revision: 3/17/98 
s UHuputs- L and Pen from the collision phase scripts 
E ots 


a ay 


Flags corresponding to the cargo compartments that will 
release oil 


= 5 
while L > IBH(j) 

A 
end * at this end, BH(j) is the first bulkhead aft of 
the damage location 


INIM Te 
ie UO EU Sea ee 
Tan ease 
Tank1BOT = 1; 
1f Pen »-24 


Tank ioa- qu 
end 
end 
end 
if F. SS) 3¢ 
ele? 9 
Tank2P =< T: 
TankZ2207 i 
if Pen >=24 
Tarso tele 
end 
end 
end 
TR NU 
Lt roptUre2? > 0 
NOUS Pom 1; 
Tank3BOT = 1; 
if Pen >=24 
Tank So ie 
end 
end 
end 
LeS coe. 


end 
a 


end 
1E] 


end 
ali 


end 


o ec 
Tank4P = 1; 
Tank4BOT = 
if Pen >=24 


ld 


Tank48 =1: 


end 
end 


== 6; 

oc lO, 
Tank >= ik 
TanksBoTr = 
if Pen 2-24 

Tank5S 

end 

end 


m EI 

LE rupture? mm0 
mec 
Wan om 
1f Pen >=24 

Tank6s 

end 

end 


== 8; 

if ruptüre2 >70 
So li ST 
TankóBOT = 
SOIT 

SlopTk 

end 

end 


d 


OS 


1 


=l: 


Jr 


-2; 


% write.m 


Ei cumcEECOMS-uN -mmasssenedeeelculbxted resu Ca for later 


analysis. 
z ae EH e 


% Output: Generated variates and values 

% Variable label Variable description 

% vil Struck ship speed (kilik vel 
y-vel] 

% V2 Striking ship speed (kt) [x- 
vel, y-vel] 

% alpha Striking ship half-entrance 
angle (degrees) 

% theta Collision angle (degrees 
relative, [rom struck ship) 

% L Location Gf Collision point sonore rae 
ship(meters from FP) 

% E Energy absorbed by structure 
% 2 Tine Tor collis onto O ar 
% de time step for time domain 
analysis 

$ R ftag Set bou eo um 
SEE AO ELIAS 

% Pen depth of penetration during 
collision 

$ Len Length- orf collisione eS 


% Date created: 11/3/97 
% Last updated: 2/16/98 
il time <= 45 


EDSR, 
Eshell(i)=Emem; 
VIU 

NONE 

bow alpha(i) = alpha; 


if phi « 90 
LOC(1)zmax(0,Lo-((Beam*cos(alpha'/2))); 
else LOC (1) =min (LBP, Lo+ ( (Beam*cos (alpha) /2))); 


end 

DOE ae ales 

ICP(i) = Lo/LBP; 

FCP(i) = L/LBP; 

Len(i) = abs((LOC(1)/LBP)-FCP(1)); 

R(1) - rüptüre, “= flag set to 2 2e shel ene ame 
coge 

CORSI EET EE 

P(i) = Pmax; 

Gen(i) =~ (ICP) Hep O. 

Minorsky) GREC 


On 


Time i) es echoes 
Guti Low ae 
Mass(i) - ml; 
end 





% howmany.m 


o 


them from the population 
numberzeroed = 0; 
for d= rn. 
E Time (QUO OL 
numberzeroed = numberzeroed +1; 
end 
end 
numberzeroed 
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s scripto Check fork ntegration error cases, 


and remove 


% calibrate.m 


% This script is used in Che Input par roo e 


process. 
The first step is to remove cases where no outflow 
occurs from the sample 

% population. Next an error function based on a 
squares" analysis is calculated. 


oe 


"leac 


oco 


Date created: 
Last revision: 


27 ey lO 
3727/96 


% Loputs: Result vectors from write.m 

% OUI: Result vectors with non-rupture cases removed 
% Value of error function tores ccu 
simulations 

% Plots of output with zero outflow cases 
removed 


Remove non-rupture cases from sample population 


o 


% First, zero out the penetration, length of damage 
and center of damage 
% elements from the cases where no hull rupture 
occurred 
A 
Bee a) eG, 
P(i) = 0; 
Len(i) = 0; 
Con = 0, 
bow_alpha(1) = 0; 
CA(1) = 0; 
TOP OR O 
end 
end 
% Now, remove those cases from the population using the 
"nonzeros" function 
Len = nonzeros (Len); 
P = nonzeros (P); 
Cen - nonzeros (Cen); 
bow alpha - nonzeros(bow alpha); 
CA = nonmzeros (Cua 
LCP = nonzeros (ICP); 
% Record the length of the resulting vectors Lor ruci: 


use in Che error function 
s = min(length(Len), length(P)); 
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Record the simulation population in each bin 


aser the number Of bans 


iS Oe IO DIE CHO Oo wally, 
node (0.0580. 05m |= 
eS yt iF; 
nba = [5:10:175]; 


Moba = (0-25-90) > 
mbeut == 10-0501. 0. 51 


% record the populations 
fore NS CBS 

Ee] hist (Len,nbl); 
os Huet cens e 
[ 
[ 
[ 
[ 


Ho wy —g 


collang, dli = NisSC(CA, mb e 
hist(bow alpha,nbba); 
heu T CP. nbl 
Lu srioummEbc roc > 


bowang,e] 
cob e 
exuti g] 


Convert the bin populations to probability density 


tuncer S 


by o> Oe Se 

1 1 (s0708). 

C e MiS a 

collang collan scu 
bowang = bowang/(s*5) 
collpt colla Ss ooo 
Se OE ln As O 


p 


Plot the resulting pdf's 
figure (2) 
eor 
colormap (white) 
hold on 
Darr pD) 
=O TOS T5 EE CE NE 
Vol24, 9645,0256 705 S67 re aor 0 Soo Omi Olle 
POLOEN verses) 
plotoubpse. deu 
axis([0 .4 0 30]) 
legend('MARPOL Standard','Calculated Distribution',1) 
ylabel ('Probability Density') 
xlabel('Transverse Penetration (Penetration/Beam) ') 
title('Penetration of Transverse Damage PDF') 


figure?) 
er 
colormap (white) 
ne Ieee 
bar la, 1) 
LO Od E 
y=(11.95,3.5,0.35,0.35,0]; 
201 


¡A yee oe 

puot: Tie M) 

axis([0 .4 0 20]) 

vlabel ('Probability Density') 

xlabel ('Length of damaged section (Length/LBP) ') 
legend('MARPOL Standard', 'Calculated Dieceri bution | 
title ('Longitudinal Length of Damage PDF') 


figure (4) 

idus 

colormap (white) 
hold on 
Dard 2e] 
A A 
YO 
¡AA 
piano, cal 
AO: IO EOS 
Xlabel ('Longitudinal center of damaged section 


(66 / La) 


vlabel ('Probability Density') 
legend ('MARPOL Standard', 'Calculated Distribution AN 
title('Longitudinal Center of Damage PDF') 


% This part is to plot the input pdf's and the generated 
histograms.. 


46)^2 


136)” 


Frios S) 

colormap (white) 

Gite 

MOS SO 
AOS OO O 
Ss Tona =>. 

SEO A Sd 


% disti(i)= (l/sqrt (2*pi*sigma 2) ) "ezr l aE 

) / (2 mese E 

% dist2(i) = 3*(1/sqrt (2*pi*sigma”™2) ) *exp((- (i- 
2)/ (2*sitemea 2)05 

% dist (1) Stews ele NF r dee EM EE 

send 

Doc 


bar(d eed aria) 

plot (x e ar] 

pucetdnmia, co SCIL. US 

axis tici 

ylabel('Probability Density!) 
xlabel('Collision Angle (degrees relative) ') 
legend('Target PDF', 'Actual Histogram', 2) 
title('Collision Angle PDF') 


figure(6) 
colormap (white) 
CIT 

dist = Q0; 

x=- [0 eoo: 


Ee 
TORA ls ol 
dist(i)= (1/sqrt (2*pi*sigma”2) ) *exp((-(i- 

cM ey (2 * sigma 2) ) > 

end 

Douce on 

bar (e, bowang) 

planes disti. =) 

plot (nbba, bowang, 'k*') 

axis tight 

ylabel('Probability Density') 

xlabel('Bow Half-entrance Angle (degrees) ') 

legend ("Target PDF", ‘Actuarenisvogram’, 1) 

title('Striking Ship Bow Half-Entrance Angle PDF') 


figure (7) 

colormap (white) 

C LE 

Dod rc 

bons. colli) 

OO ODO e 

dust e (1,1701 

pct a E) 

ES RES MODE) 

plot (nbl coto kt 

ylabel('Probability Density') 
xlabel('Initial Collision Point (Loc/LBP)') 
legend('Target PDF','Actual Histogram',1) 
title ("Initial Pornt. ot Contact ln Colla ion men 


figure (8) 

colormap (white) 

clr 

hold on 

bar (g,outflw) 

cals ( (0 O.5. G50) 

Plot (about out iiw a e 
ylabel('Probability Density') 
xlabel l c pecil ic outtilem or] rgo 
legend('Generated Histogram',1) 
title (TOIL Out lova PDE 


% Calculate tne error functions 

BRP jeep (i) =14 98) 2 12902 75) oa MEE ee 
+ (p(4) GO. 56) 202 p(s) "2 4: p(6)° 2 42017) 2 ap eee 
9) eo eon hover 

EFP = EFP * p(11)^2 «p(12)^2 +p(13)"2 sp(i4) 2 +p(15) "2 
+p(16) 2 TT o omo) co c DU DE 

EFP = I- sor E 20: 

A O A E A SAS a ee 
DIS ee A 1375). A O Ns A a 
4 AAA A 008509.) 0 A os) onan s 


Duo IS) S : + 1(15)°2 + 1(16) 2X1 02 OE MEM 
A C20) 
EFL = 1- ld O 
EFC = (c(1)-1) 2 + (CM) 112 FI A 
1)^2 + (e(5) = 1)°2 + (c6) 1) 25 005 ee 
a L 


sop or cO D MES 2 + (c(10) =- 1) 2; 


EFC = 1-sqrt (EFC)/10; 


AVE = (EFP+EFL+EFC) /3; 
Product = EFP*EFL; 


EFP, EFC, EFL, AVE, Product 


Pzero — "(Ol 
fOr j= eee 
olen o = — 


* Calculate probability of zero outflow, and mean outflow 


Pzero = Pzero +1; 
end 
end 
Pzero=Pzero/s 
Mean = sum(Outflow)/s 


% output.m 


oe 


A script to generate all output graphics from a run 


Date created EN 25 S 


o? oe 


Last revision: 3/27/98 
% eu Result vectors from write.m 
% uU Plots of input and result pdf's for this set 
of simulations 
% Record the length of the resulting vectors for future 
use 


S - min(length(P), length(Len)); 


% Record the simulation popúlation in cachi pin 


5$ set the number of bins 


o 


11599857 500:05.: 02:05 $1] 4 * number of bins for 
penetration 

nol = 002S 0.05 T975]. %7 NUME ROL DINER o 
length 

De eee OO se UE v number ot balms tor Nc Mer OE 
damage 

nba te 5e T7599 t5 number of bins for collision 


angle 
nbba = [2.5:5:87.5] ; 
Hooue = 1020 20.045 e 
outflow 


number of bins for bow angle 
number of bins for vous 


o? o 


% record the populations 

ee hist (P,nbp) ; 

ae hist (Len, nbl); 

[c,z] = hist (Cen,nbc) ; 

(collang,d] =. hist (CA, nbad); 
[bowang, e] hist (bow alpha,nbba) ; 
(Cola ter ste LC Abi 

[out El wa) Nist (Out low, neeourt e 


% Convert the bin populations to probability density 
functions 


p / (80.05) ; 

ii (SOS E 

c =a cb) 

COUP Nc ENSE 


p 
l 


Ho gu 


bowang = bowang/(s*5) ; 

cojus c MEA cOMUE e 

out Sodi Soi 
% Plot the resulting pdf's 
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figure (2) 
(ud: 
colormap (white) 
hele yon 
bar (x,p) 
x= [07 054. Ss: o MENU ME 
Y=|(24.96, 5,07 236.0 3567,0 7567 0250 tres > One 
plotio va ae, 
plot(nbo) DOS SEE 
axes GN S EP SON 
legend('MARPOL Standard','Calculated Distribution!) 
vlabel ('Probability Density') 
xlabel ('Transverse Penetration (Penetration/Beam) ') 


figure (3) 

QE 

colormap (white) 

holla on 

bar (oye 

Se LO eee 3. E 

v= (21295325002), Cfo Ole 

plc c 2M) 

DIOR abit hea) 

axist[ow 5 0M NN 

yláabel( l Probabili EDENS IEN 

xlabel('Length of damaged section (Length/LBP) ') 
legend('MARPOL Standard', 'Calculated Distribution',1) 


figure (4) 

ee 

colormap (white) 

Hod orn 

lam m) 

= Oy 1 EBD 

y=[1,1,0]; 

DUOe toe ce, es) 

PLOCA Et 

axioo T O20 oss) 

axis tight 

xlabel ('Longitudinal center of damaged section 
(boc DEB 

ylabel('Probability Density!) 

legend('Calculated Distribution','MARPOL Standard',1) 


* This part is to plot the input pdf's and the generaten 
histogr arona 


figure) 

colormap (white) 

Ale 

scs oS ed 

dist = D9O0U5S. 0.0056 5 
SLOPE ES 


206 


% Geel (Jana / sort (2*o1 “sigma ate eo 
46)^2)/(2*sigma^2)) 

x gus (al). 
13569 2)/(2*sigma 2) 


= 3* (1/sqrt(2*pi*sigma^2))*exp((-(i- 
S 
% Gast a 


Dard Collang) 

pec cisco ny 

Pier (nba, eellang, "k= 

axis tight 

ylabel('Probability Density') 
xlabel('Collision Angle (degrees relative) ') 
legend ('Target PDF', 'Actual Histogram' , 2) 


figure (6) 
colormap (white) 


gus 

dist = 0; 
x= DEIS 
sigma=5; 


ju TE E 
dist (i)=.(1/sart(2*pi*sigma 2))ʻexp( (1 

39)^2)/(2*siqma^2)); 

end 

ones 

bar (e, bowang) 

PLOT Si toe) 

plot (nbba,bowang, 'k*') 

axis tight 

ylabel('Probability Density') 

xlabel('Bow Half-entrance Angle (degrees) ') 

legend('Target PDF', 'Actual Histogram',1) 


figure (7) 

colorma” (white) 

Cisne 

Nolde n 

pori col 

so ore QU Qs 

dise A UM 

pios M NM) 

axis POE O o 

a e 

Pim le OLMO CA an 

vlabel ('Probability Density') 

xlabel(' initial Collision Pointe. (hoc, ibm. 
legend ('Target PDF','Actual Histogram',0) 


figure (8) 
colormap (white) 
pq 
hold on 
bar(g,outflw) 
207 


taxis ()Omemse 0 Siig 
¡A NONE. 

axis([sS9MOo MESE NCMO 

axis tigne 

ylabel('Probability Density') 

xlabel ("Specific Outflow (ontario 
legend('Generated Probability Density',1) 
titled] OU Low ERE 


* Calculate probability of zero outflow, and mean outflow 
Pzero = 0; 


for JEn 
ice Oou C laa a 
Pzero = Pzero +1; 
end 
end 
Pzero=Pzero/s 
Mean = sum(Outflow) /s 
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%  joint.m 


% Ars crio co lgenerate JOlMs probability tuner 
Gea ST Rr orm 


Date Created 2 25 9S 


o? o? 


Last- r- vicion. 3/27/98 
$ puts: Result vectors Eren write i 
% QUE SUE Plots Joint pdf's for this set on smile ome 
Z =) Zeros 20 20I * Sets up space for the joint pdf matrix 
values 


toral Ot pehem lenin Ane Senekiza 6 tome paises 
jog I cec 


$ Start at the first bin, and find the appropriate length 
and penetration bins 
IN lier 
I 
while Len(i) » nbl(k) 
| M 
end 
while P(i) > nbp(m) 
meum de 
end 


% Add one to the current bin population 
A ( (0 99) eh eR 
end 


Reset bin values to produce a "Sensible" plot 
orem qo 
nbp (m) 

mod me 
d 


|. dg» 

e 
0 0O 
aii 
m 
I I 

oS 
Se 
Ul Ul 


o? o9 o? o9 o9 


D 


= 


cor expo cC 
22/39 9050705), 


N o 


o 


+ and plot 

figure (9) 

colormap (jet) 

alr 

MeshoGisie (mite nao a. 

mesh (nbl,nbp, Z) 

View (an 

axis (iG. CD EDEN 

axis tight 

ylabel ('Longitudinal Extent') 
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Zlabel('Probability Density') 


'Transverse Extent') | A 
ic PDF for Transverse and Longitudinal Extents of 


Damage') | 
colormap (lines) 
brighten(-0.3) 
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